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The late-stage demixing following spinodal decomposition of a three-dimensional sym- 
metric binary fluid mixture is studied numerically, using a thermodynamicaly consistent 
lattice Boltzmann method. We combine results from simulations with different numerical 
parameters to obtain a unprecendented range of length and time scales when expressed 
in reduced physical units. (These are the length and time units derived from fluid density, 
viscosity, and interfacial tension.) Using eight large (256 3 ) runs, the resulting composite 
graph of reduced domain size I against reduced time t covers 1 < ^ < 10 5 , 10 < t < 10 s . 
Our data is consistent with the dynamical scaling hypothesis, that l(t) is a universal 
scaling curve. We give the first detailed statistical analysis of fluid motion, rather than 
just domain evolution, in simulations of this kind, and introduce scaling plots for sev- 
eral quantities derived from the fluid velocity and velocity gradient fields. Using the 
conventional definition of Reynolds number for this problem, Re^ — l&l/dt, we attain 
values approaching 350. At Re^ > 100 (which requires t > 10 6 ) we find clear evidence 
of Furukawa's inertial scaling (I ~ t 2 ^ 3 ), although the crossover from the viscous regime 
(I ~ t) is both broad and late (10 2 < t < 10 6 ). Though it cannot be ruled out, we find 
no indication that Re^ is self-limiting (I ~ t 1 / 2 ) at late times, as recently proposed by 
Grant and Elder. Detailed study of the velocity fields confirm that, for our most inertial 
runs, the RMS ratio of nonlinear to viscous terms in the Navier Stokes equation, R2, is 
of order ten, with the fluid mixture showing incipient turbulent characteristics. However, 
we cannot go far enough into the inertial regime to obtain a clear length separation of 
domain size, Taylor microscale, and Kolmogorov scale, as would be needed to test a re- 
cent 'extended' scaling theory of Kendon (in which R2 is self- limiting but Re$ not). To 
obtain our results has required careful steering of several numerical control parameters so 
as to maintain adequate algorithmic stability, efficiency and isotropy, while eliminating 
unwanted residual diffusion. (We argue that the latter affects some studies in the litera- 
ture which report I ~ t 2 ! 3 for t < 10 4 .) We analyse the various sources of error and find 
them just within acceptable levels (a few percent each) in most of our datasets. To bring 
these under significantly better control, or to go much further into the inertial regime, 
would require much larger computational resources and/or a breakthrough in algorithm 
design. 

f Present address: Optics Section, The Blackett Laboratory, Imperial College, London, SW7 
2BW, UK. 
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1. Introduction 

Spinodal decomposition occurs when a fluid mixture of two species A and B, forming 
a single homogeneous phase at high temperature T, undergoes spontaneous demixing 
following a sudden drop in temperature (or 'quench'). For suitable compositions and 
quenches, one enters the 'spinodal' regime in which the initial homogeneous phase is 
locally unstable to small fluctuations. (Elsewhere one finds instead a nucleation and 
growth mechanism which is not the subject of this paper.) For compositions close to 
50/50, there then arises, after an early period of interdiffusion, a bicontinuous domain 
structure in which patches of A-rich and i?-rich fluid are separated by sharply defined 
interfaces. The sharpness depends on the temperature drop; we assume a 'deep quench' 
for which the interfacial thickness is, in practice, on a molecular rather than macroscopic 
scale. In this late-stage structure, the local compositions of the fluid patches correspond 
to those of the two bulk phases in coexistence; the interfacial tension approaches <r, 
its equilibrium value. Although locally close to equilibrium everywhere, the structure 
then continues to evolve so as to reduce its interfacial area. Local interfacial curvature 
causes stresses (equivalently, Laplace pressures) to arise, which drive fluid motion. The 
interface then evolves smoothly with time between isolated 'pinchoff events' or topological 
reconnections. In principle these events reintroduce molecular physics at the short scale; 
however it is generally assumed that pinchoff, once initiated, occurs rapidly enough not to 



impede the coarsening process (but see Jury et al. 19995; Brenner et al. 1997). Likewise 



it is usually assumed that at late times, the presence of thermal noise in the system is 
irrelevant, at least for deep quenches (but see G Gonnella fc Ycomans 1999] ): the problem 



is thus one of deterministic, isothermal fluid motion coupled to a moving interface. Precise 
details of the random initial condition, which is inherited from the earlier diffusive stage, 
are also thought to be unimportant (assuming that no long-range correlations are initially 
present). 

For simplicity we address in this paper only the maximally symmetric case of two 
incompressible fluids with identical physical properties (shear viscosity ry, density p), and 
also equal volume fractions, that have undergone a deep quench. With the assumptions 
made above, all such fluid mixtures should, in the late stages, behave in a similar manner. 
More precisely, the dynamical scaling hypothesis is that, if one defines units of length and 
of time by 

L = V 2 /(pa) , T = n 3 /(po- 2 ) (1.1) 

(which are the only such units derivable from 77, p, a) then at late times any characteristic 
structural length L(T) should evolve with time T according to 

dL/ dT = (Lq/Tq) ip(L/ Lq) (1.2) 

where <p{x) is the same function for all such fluids. (A specific choice of definition for L is 
made later on, in terms of the mean domain size.) Integrating this once gives a universal 
late-stage scaling 

l = l{t) (1.3) 
where we introduce 'reduced physical units', 

l = L/L , t= (T — T mt )/T . (1.4) 

Here Ti nt an offset that is nonuniversal: it depends on the initial condition as fixed by 
the early stage diffusion processes. (Note that in this paper, the symbol t is reserved for 
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the reduced physical time; unsealed time is denoted T, and temperature T. An overdot 
means time derivative in whatever units are being used.) 

The form of l(t) has been discussed by several authors, notably Siggia (1979) and 
Furukawa (1985| ). Siggia argued that, for t <C 1, the interfacial forces induce a creeping 



flow of the fluid; simple force balance in the Navier-Stokes equation then gives I oc 
t in this, the 'viscous hydrodynamic' regime. Note that, in a creeping flow, the fluid 
velocity depends only on the instantaneous structure of the interface. (This is why the 
nonuniversal offset is in T and not L.) At later times, the force balance was argued by 
Furukawa to entail viscous and inertial effects; balancing these giv es I oc t 2 ^ for t ^S> 1, 
the 'inertial hydrodynamic' regime. It has recently been shown by Kendon (2000 ) that 
Furukawa's assumption of a single characteristic length (for velocity gradients as well as 
interfacial structure) is inconsistent with energy conservation; her more detailed analysis 
nonetheless recovers I oc i 2 / 3 for the domain size. Kendon's arguments, with those of 
Siggia and Furukawa, are discussed in § ||, [|. 

For a general review of late-st age spinodal decomposition and other aspects of phase 
separation kinetics, see that of Bray (1994). The problem is clearly intractable ana- 
lytically: it involves a moving boundary with a complicated and non-constant topology 
whose initial condition is defined, implicitly, by the preceding, early-time diffusion. These 
features render it equally intractable to many numerical algorithms that might perform 
well for other fluid mechanics problems. Indeed, symmetrical spinodal decomposition has 
become a benchmark for various so-called 'mesoscale' simulation techniques, developed 
to address the statistical dynamics of fluids with microstructure. The results from dif- 
ferent techniques can be compared, not only with each other and with experiment (with 
the caveat that one cannot realize exact symmetry between fluids in the laboratory), but 
with the predictions of the various scaling theories already mentioned. 

In the present work, we study in detail the physics of spinodal decomposition for a sym- 
metrical binary fluid using the Lattice Boltzmann (LB) technique (Higuera et al. 1989), 
in a therm odynamically consistent form pioneered by the group of Yeomans (see Swift 
et al. 1996 ). Our work, of which a preliminary report appeared in Kendon et al. (1999 ), 
advances significantly the state of the art for simulations of spinodal decomposition, and 
for LB simulations of fluid mixtures. In any such simulation, a balance must be struck 
between discretization error at small scales, and finite size errors (arising in our case 
from periodic boundary conditions) at large ones; this compromise is quite subtle, as we 
discuss below. It means that any individual simulation run can produce only around one 
decade of data for the l(t) curve. This is true for all first-principles simulation methods: 
in three dimensions there cannot be more than two decades, or at most three, sepa- 
rating the discretization length from the system size, before deduction of a half-decade 
safety margin at each end. (Three decades before such deduction is optimistic; it means 
simulating at least 10 9 degrees of freedom.) 

Despite this restriction, by careful scaling and combination of separate datasets for 
eight large (256 3 ) simulation runs, we are able to access an unprecedented range of I 
and t (five and seven decades respectively) including regions of the 1(f) curve not studied 
previously. We find nothing to contradict the universality of Equation ( |l.3| ), but nor 
can we completely rule out violations of it. We gain the first unambiguous evidence for 
a regime in which inertial effects dominate over viscous ones, and find clear evidence 
for i 2 / 3 scaling in this regime. In this and other regions of the 1(f) curve, we study 
the statistics, not only of the interfacial structure, but also of the fluid velocity. (The 
latter was not addressed in detail by previous simulations.) On entering the inertial 
hydrodynamic regime we find some evidence for breakdown of simple scaling of velocity 
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gradients, as proposed by Kendon (200C), but our data does not extend far enough into 
this regime to offer a meaningful test of her alternative proposals. 

To obtain our new results, we have had to push the LB technique to its limits. For 
statistics based on the domain size, errors at the level of several percent, arising from each 
of several different sources (residual diffusion, lattice anisotropy etc.) remain. We do make 
a systematic attempt to identify and minimize the various sources of errors - a somewhat 
arduous task that, our work suggests, has been neglected in several previous studies. 
The errors for some of our velocity statistics (especially those for spatial derivatives of 
the velocity) are much larger. Nonetheless we present the data, such as it is, because 
it highlights several issues both in the physics of spinodal decomposition and in how 
simulation results should be obtained, analysed and interpreted. 

The rest of this paper is organized as follows. § ^ outlines the thermodynamics of the 
binary fluid system, and § |^ its governing equations. § |4| and § |^ outline the simple and 
extended scaling analyses referred to above. § ^describes the LB method in the form that 
we use; § [?] describes how the simulation parameters are chosen. § || outlines a number 
of validation tests. § ^ gives our results for the evolution of the interfacial structure, 
[To| those for the velocity field and § [□] those for the velocity derivatives and related 



quantities. § 12 summarizes our conclusions. Two appendices give further information on 



the effects of residual fluid compressibility in the LB method and on the relation between 
our work and that of previous authors. 



2. Thermodynamics 

Although we are interested in the late-stage demixing of two isothermal, incompress- 
ible fluids separated by sharp interfaces, the LB method resorts to a more fundamental 
approach, in which these interfaces are described as excitations of a thermodynamic field 
theory. The central object is the Helmholtz free energy 

F = E-TS 1 (2.1) 

where E is the internal energy, T the temperature and S the entropy of the system. 

In a system at fixed volume V , and fixed contents and temperature, equilibrium states 
are given by global minima of the free energy, F . For a symmetric fluid mixture, F is a 
functional of a single composition variable <j>(r), defined as <j) — (ua — «s)/(^A +hb) 
where the n's are number densities, and of the mean fluid density p = ua + tib- (We 
take unit mass for A and B particles without loss of generality.) In the incompressible 
case, p is fixed; we leave it as a parameter in what follows. Further restricting attention 
to homogenous states (so that (j) is the same everywhere), we can write 

F/V = V(0). (2.2) 

Within mean-field theories of fluid demixing, one predicts that V has everywhere positive 
curvature at high temperatures, but becomes concave below a critical temperature T c . 
The resulting curve is as shown in Figure [l], with symmetric minima at ±(f>*. Below T c , 
the free energy is therefore minimized by creating two bulk domains (of equal volume) 
at compositions ±cf>* instead of a single homogenous phase with </> = (which is our 
presumed initial condition). The same phase separation occurs for any other <f> between 
±0*, but in this case the domain volumes are unequal; for sufficient asymmetry this 
causes depercolation. (In a depercolated, droplet structure, coarsening can only occur by 
diffusion so that the scaling arguments given above cease to apply. We do not address 
this here.) 

The resulting phase diagram is shown in Figure 0. Spinodal decomposition occurs for 
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Figure 1. Model potential for phase separation, a symmetric double well, V(4>). The 
equilibrium values of the order parameter are ±0*. 



QUENCH 



mixed binodal line 

spinodal line 




Figure 2. Phase diagram for spinodal decomposition. The order parameter is 
<f> = {n.A — n_g)/(n.4 + ub) with p — ua + ub the mean fluid density. The temperature 
axis shows the critical temperature, T c , below which the system starts to separate, and above 
which it remains completely mixed. 



any quench that leaves the system beneath the spinodal line, on which d 2 V/d<fi 2 changes 
sign. Immediately after such a quench, the system is locally unstable: the free energy 
can be lowered, in any local neighbourhood, by creating two domains whose composition 
differs only infinitesimally from the initial one. (The resulting free energy density lies on 
a line connecting two points on V(4>) at the new compositions; in the convex region, this 
causes a reduction in F.) Accordingly, infinitesimal fluctuations will grow by diffusion 
until there is local coexistence of domains at compositions approaching ±<f>* . 

To describe quantitatively both the domains and the interfaces between them, one 
must specify not just V(<fi) but the free energy functional, F [</>]. An acceptable choice is 



the square gradient model (see Bray 1994) 



F[<f>] = |dr {V(0) + ||V0| 2 }, (2.3) 

where V(<f>) is as shown in Figure and the term in k penalizes sharp gradients in 
composition. This ensures smooth local deviations from ±0* near the interface, and 
provides a nonzero interfacial tension a which can be calculated as follows. We consider a 
flat interface between two domains, introducing a coordinate normal to it, g. Stationarity 
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of F requires 



HV 2 4> — K—^r = 



d 2 4> _ dV(^) 



(2.4) 



Integrating this once across the interface and setting g — 0, <j) — at its centre gives 



The excess free energy per unit area of interface is then given by 



a = / dg 



k ( d(p 



(2.5) 



(2.6) 



2 \dg, 

By exploiting the fact that — * in the bulk fluid, and using Equation (2.5), we obtain 

d0(2 K ) 1/2 [V(0)-V(r)] 1/2 . (2.7) 



Given a form for the potential, V((f>), a value for the interfacial tension can thus be 
calculated. This is done for the model used in our simulations in § ^. 

We now turn to the (exchange) chemical potential, fj,, which describes the change in 
F for a small local change in composition: 



SF dV „o 



(2.8) 



Within the LB approach, the coupling between interfacial and fluid motion arises as 
follows. In the presence of a nonuniform composition, there is a thermodynamic force 
density — ^>V/i acting at each point on the fluid. (The two species are pulled in opposite 
directions by the chemical potential gradient; the net force vanishes only if <f) = 0.) This 
force density can also be written as the divergence of a 'chemical' pressure tensor: 

</>V/i = v.v chem 



where it is a straightforward exercise to confirm that 



'a/3 



Sat 



dV 



V- K {0V 2 0+i|V^| 2 } 



K(Va0)(V/J0). 



(2.9) 



(2.10) 



Note that only the last term is anisotropic; the rest contributes, in effect, to the isotropic 



fluid pressure P. By integrating (2.9) across an interface, and using Eq. (2.1C), one finds 
that there is, in static equilibrium, a finite pressure difference across a curved interface, 
called the Laplace pressure: 

AP = a/C (2.11) 

where K, is the interfacial curvature. 

Throughout the above, our description in terms of a smooth composition variable 0(r), 
usually known as the order parameter, assumes a coarse-graining so that the smallest 
length scale under consideration is larger than the average distance between molecules. 
In equilibrium, this coarse-graining is an almost trivial operation, but for the dynamical 
description desired below, certain conditions must be met. On scales smaller than the 
coarse-graining length, the system must remain in local equilibrium, while the variations 
of interest at larger scales must be slow on the scale of the time it takes for that local 
equilibrium to be reached. 

This does not mean that the microscopic scales can be forgotten from here on. Although 



Inertial effects in spinodal decomposition 



7 



usually a macroscopic description is sufficient to fully understand the system, ultimately 
it is still the microscopic interactions that are driving the system and determining the 
dynamics. In particular, any interface between the two fluids will have a microscopic size 
and structure. It is always a possibility that the microscopic behaviour can intrude at 
the macroscopic level (for example, by interfering with pinchoff) and change the results 
predicted by any simple macroscopic considerations. In particular, when using numerical 
models, care must be taken that the microscopic behaviour in these models is admissible. 



3. Governing Equations 

The equation of motion for <p is taken to be a convective diffusion equation of Calm 
Hilliard type (see [Bray 1994 |Swift et al. 1996| ) 

4> + v.V0 = AfV 2 At = -A/V 2 | K V 2 - , (3.1) 

where M is an order-parameter mobility (here assumed independent of <f>) that controls 
the strength of the diffusion, and v(r) is the fluid velocity. This equation states that the 
order parameter responds to composition gradients by diffusion (the MW 2 fi term), and 
also changes with time because it is advected by the fluid flow (the v.V(/) term). 

The fluid velocity in turn obeys the Navier-Stokes equation (NSE), which for an in- 
compressible fluid reads 

p[v + (v.V)v] =r/V 2 v- V.V th . (3.2) 

Here V^g is the 'thermodynamic' (or conservative) part of the pressure tensor, and con- 
tains two pieces: an isotropic contribution P5 a p, chosen t o ma intain constant p, an d the 
'chemical' pressure tensor, V^ em , defined previously in ( |2.10 ). Recall that by fl2.9|) , the 
chemical term —^/p chem can equally well be represented as a body force density — </>V/i 
acting on the fluid, so that Eq. (3.2) also reads 

p [v + (v. V)v] = r]V 2 v ~ 4>Vn - VP . (3.3) 

Within the LB approach of Swift et al. (1996| ), the governing equations are solved by 
relaxing slightly the requirement of fluid incompressibility. We return to this in § |[ 



4. Simple Scaling Analyses 

The pair of coupled nonlinear differential equations, (|3.l|) and (3.2), are intractable, 
but various dimensional and scaling ideas may be used to find out how fast the domains 
grow once the diffusive period is over. All these analyses assume that the interface can 
be characterized by a single length scale — that is, it is basically smooth, with radii of 
curvature that scale as the domain size itself, which is much larger than the interfacial 
thickness. 

Many domain-scale length measures are possible; we use L(T), the inverse of the first 
moment of the spherically averaged structure factor, S(k,T): 

L(T) = 2 "/fcS(fc,T)dfc' 
where k = |k| is the modulus of the wave vector in Fourier space, and 

S(k,T) = (cj>(k,T)<p(-k,T)) (4.2) 
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with 0(k, T) the spatial Fourier transform of the order parameter. The angle brackets 
denote an average over a shell in k space at fixed k. 

The aim of scaling analyses is to find the form of the time dependence of L(T) by 
considering the NSE (3.3), and balancing the force from the interface,— <f>Vfi, against the 
viscous and inertial terms which tend to oppose its motion. The interfacial force density, 
— 0V/i, can be approximated as follows. The curvature, JC, is of order l/L, since L(T) 
is roughly the size of the domains. This sets the scale of •p chem through (2.11), as a/L. 
Likewise the gradient operator, V, can be approximated by 1/L(T) in Equation (2.10), 
which then reads 

-0V M ^. (4.3) 



Now we turn to the remaining terms in the NSE (3.3). We start by assuming that the 
length L also controls the magnitude of V as far as velocity gradients are concerned. 
Approximating also the fluid velocity, v, by the velocity of the interface L(T), gives for 
the viscous and inertial terms respectively 



ryV v ~ r\ 



L 

I 2 ' 



£2 

p[v+(v.V)v] ~ P L + p— . 



(4.4) 
(4.5) 



Under conditions in which the inertial terms are negligible, the force from the interface 
will be balanced by the viscous force, giving L/L 2 ~ a/(r]L 2 ). Integrating this gives, 



L ~ - (T - T mt ) 



(4.6) 



Thus the domain size is predicted to grow linearly with time in t he region where the 
fluid flow is viscous dominated. This is the r esult of [Biggia (1979| ). Linear growth ha s 
been reported in experi ments by, for example, Kubota et al. (1992 ); Chen et al. (1993 ) 
Hashimoto et al. (1994), and in simula ti ons incorporating hydr o dynamics by Koga fe 



Bastca fe Lebowitz (1997| ) and |Jury et al. (19996 ) 



Kawasaki (199l|);|Puri fc Dunweg (1992|); [Alexander et al. (1993[ ); |Laradji et al. (1996| ); 



To find the growth rate in the inertial region, Furukawa (1985 ) balanced instead the 
inertial and interfacial terms; assuming again only one relevant length, he obtained 



T 



a 1 



(4.7) 



Integrating this twice gives, for large enough T, L 3 ~ aT 2 /p, so that the domain size 
grows as L ~ T 2 / 3 . This result has not yet been observed experimentally (for reasons we 
discuss later, § 12). There are a few previous claims to see this in simulation (Ma et al. 
Appert et al. 1995; Lookman et al. 1996), but none reliably establish dominance 



1992 



of inertial over viscous forces as we do below in 



10 



Comparing the results of (4.6) and (4.7) allows us to estimate a characteristic domain 
size, L = L* , and characteristic time, T = T* + Ti nt at which the crossover from viscous 
to inertial scaling occurs. (To be precise, we can define L*,T* by the interception of 
asymptotes on a log- log plot.) This leads to L* ~ Lq, T* ~ To, with Lo, Tq defined in 1.1. 
Converting to reduced physical units / = L/ Lq, t = (T — Ti nt )/To as defined previously, 
and invoking the dynamical scaling hypothesis, we have respectively 



I = bit 



I 



52/3 



2/3 



t<f viscous regime 
t > t* inertial regime, 



(4- 
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where 61,62/3 and t* are dimensionless numbers that should be universal to all incom- 
pressible, fully symmetric, deep quenched fluid mixtures. What scaling theories cannot 
predict, of course, are values of the universal constants 61,62/3,4*, other than to state 
that these are 'of order unity'. 

In fact our simulations show that t* = T*/Tq is between 10 4 and 10 5 , which is 'of 
order unity' only in a rather unhelpful sense; the implications of this are discussed in § 
|l2| below. We will also find that the crossover region, between the asymptotes described 
by (4.S), is several decades wide. Although there is no explicit scaling prediction for the 
behaviour within this crossover region, its width means that each individual simulation 
dataset, either within or outside the crossover, can be well-described by a single scaling 
exponent, a, such that 

I ~ b a t a , (4.9) 

where 61 ^ b a $J 62/3, and 1 ^ a ^ 2/3. We use this form below, when analysing our 
numerical data. 



5. Extended Scaling Analysis 

In what follows we will find it useful to compare directly the relative magnitude of 
the various terms in the NSE. Two ratios have therefore been defined, the RMS ratio i?i 
between the acceleration term and the viscous term, 

R 2 = OH 2 ) f5 n 

and the RMS ratio R2 between the nonlinear term and the viscous term, 

^ = ^V^vl 2 ) } ■ (5 - 2) 

Here () denote spatial averages (though an ensemble average might be preferable under 
some conditions) . These ratios obey R\ , i?2 ^> 1 where the inertial terms dominate 
and i?i,i? 2 <C 1 where the viscous term dominates. The ratio R2 we identify as the 
'true' Reynolds number, that is, a dimensionless measure of the relative importance of 
nonlinearity in the NSE. 

When i? 2 , ( |5.2D , is simplified using the normal scaling assumption (V ~ l/£), the 
result is the following estimate: 

r 2 ~ P y2 / L ~ (5.3) 
r/v/L 2 rj/p 



Assuming that the characteristic velocity scale is v ~ L, one finds (following Furukawa 



1985; Grant & Elder 1999) the Reynolds number estimate Re q 



R 2 ~ R e<j) = ^L{T)L = II = ab 2 a t 2a - 1 . (5.4) 

This 'order-parameter Reynolds number' has the advantage, in simulations, that it is 
computable from L(T) without direct access to any fluid velocity statistics. However, 
Re<f, is only a good estimate of the true Reynolds number i?2, if the simple scaling for 
velocity gradients of V — > 1/L does indeed hold. This assumption leads to the following 



paradox, as noted by Grant & Elder (1999). According to (5.4), if in the inertial region 
a = 2/3, then Re^ ~ i ia ~ x which becomes indefinitely large as time proceeds. Grant and 
Elder argued that this could not be physical, on the grounds that an infinite Reynolds 
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number would imply turbulent remixing of domains: this would limit the domain growth 



to such a speed that Re^ remained bounded at late times. (5.4) then demands a lower 
asymptotic growth exponent, a ^ 5, as t — > 00. 

However, a closer look at the scaling of all the terms in the NSE admits an alternative 



resolution. Kendon (2000) pointed out that a minimally acceptable scaling theory should 
allow, not only force balance in NSE, but a balance of terms in the global energy equation, 
which for an isothermal, incompressible fluid reads 

^ = + (5.5) 

where e = rj(ft/ a vp)(V aVp)} is the dissipation rate, and £i n is the rate of energy transfer 
from the interface to the fluid. Retaining the assumption that the interface (as opposed 



to the velocity field) has just one characteristic length, e m is readily estimated from (4.3) 
as at/L 2 . 



Applying the simple scaling for velocity gradients (V ~ l/£) to each term in (5.5) 
gives the following energy 'balance' in the inertial regime where L ~ T 2//3 : 

- pT~ 5/3 ~ -rjT- 2 + aT~ 5/3 , (5.6) 

where factors of p, T), a are retained to aid identification of the terms. At first sight 
this suggests a balance of interfacial and inertial terms, with the viscous contribution 
negligible, at late times: this is Furukawa's assumption. However, the signs show this 
to be inconsistent: the kinetic energy and the energy stored in the interface are both 
decreasing with time, so these cannot properly be balanced against each other. 

This exposes a central defect of the simple scaling analysis in the inertial regime. It is 
well known, of course, that even the simplest theories of fluid turbulence entail several 
length scales (whereas more modern 'multiscaling' theories have, in effect, infinitely many, 



Frisch 1995). In the simplest, Kolmogorov-type approach (see Frisch 1995; Kolmogorov 
1941), the important lengths are the Taylor microscalcjj] 

A = (5f](v 2 }/e) 1/2 , (5.7) 

characteristic of velocity gradients, and the Kolmogorov (dissipation) microscale, 

X d ee 2^ 3 /A) 1/4 , (5.8) 

the length scale below which nothing interesting occurs. (Energy is dissipated at or above 
the scale Arf.) 

Kendon (2000) argued that in a binary fluid system where the fluid motion has become 
turbulent, the velocity follows the interface and scales as L, but the first and second 
gradients of velocity have scalings set by A and A^ respectively, rather than by L. Within 
this simplified (Kolmogorov-level) description, the only scalings for these three lengths 
found physically admissible by Kendon for the inertial regime were A ~ T 1 ' 2 ,^ ~ 
T 5 / 12 ,L - T 2 / 3 . In the NSE, this gives for the acceleration, convection, viscous and 
driving terms the following scalings: 

pT -4/3 + pT -7/Z „ vT -7/6 + aT -i/3 (5 Q) 

The predicted outcome is thus a balance between the nonlinear and dissipative forces 
that is decoupled from the interfacial motion, while interfacial stresses balance fluid 
acceleration. The existence of a nonlinear/ viscous balance implies an asymptotically finite 

f The prefactors in the definitions of the Taylor and Kolmogorov microscales differ between 
sources. 
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Quantity 


viscous 


inertial rej 


pon 




region 


si m nip 


new 






scaling 


scaling 


L(T) 


1 


2/3 


2/3 


X = (5n(v 2 )/ £ y/ 2 


1 


2/3 


1/2 


x d = 2^ 3 /A) 1/4 


1/2 


1/2 


5/12 


V 





-1/3 


"1/3 


pv 


= 


-4/3 


-4/3 


pv.Vv 


= 0, -1 


-4/3 


-7/6 


r/V 2 v 


-2 


-5/3 


-7/6 


4>Vn 


-2 


-4/3 


-4/3 


Re^, = 11 


1 


1/3 


1/3 


Ri 


= 


1/3 


-1/6 


R 2 


= 


1/3 





e = 77(Vv) 2 


-2 


-2 


-5/3 



Table 1 . Sum mary of predic ted scaling exponents for the viscous and inertial regions. The 'new' 
scaling theory (Kendon 2000) gives the same predictions as simple scaling for the viscous region, 
except for the JNSE term pv.Vv. Entries are powers of T, an entry of indicates the quantity 
is constant, while an entry of = indicates the quantity is assumed to be zero in the viscous 
approximation . Bold entries indicate scaling predictions that differ from the simple theory. From 
Kendon (20o"o|). 



value for the ratio of the corresponding terms in NSE, that is, a finite asymptote for the 
true Reynolds number i?2- On the other hand, since the result for the domain scale, L ~ 
T 2 / 3 , survives unaltered from the simple scaling theory, the Reynolds number estimated 
from the order parameter, Re,/,, continues to grow indefinitely. This suggestion, although 
speculative, appears to resolve the issue raised by Grant and Elder, without requiring a 
change in the domain scale growth law (nor any breakdown of universality of (1.3)). 

To summarise, Kendon (200C) predicts a balance in which energy is first transferred 
from the interface (— 0V/i) to large scale fluid motion (pv). The nonlinear term (pv.Vv) 
then transfers the energy from the large scales down to smaller scales where dissipative 
forces (?7V 2 v) finally remove it from the system. The resulting energy cascade thereby 
decouples the energy input scales from the dissipation scales — a familiar enough idea in 
turbulence theory (Kolmogorov 1941). In contrast, in the viscous hydrodynamic regime, 
the simple (one-length) scaling theory is already consistent with energy conservation, and 
all its results are recovered. Kendon's predictions for scaling are summarised for ease of 
reference in table [I]. 



6. Numerical method 



The model system described by (3.1) and ( |3.2| ) was simulated numerically using a mod- 
ular LB code called Ludwig, described in detail in Desplat, Pagonabarraga & Bladon 
(2000a). It has both serial and parallel versions; the parallel code uses domain decompo- 



sition and the MPI (message passing interface) platform. Any cubic lattice can be used 
with the Ludwig code; the lattice parameter is taken as unity, as is the timestep AT, 
thereby defining 'simulation units' of length and time. 

Here we chose the D3Q15 lattice, a simple cubic arrangement in which each site com- 
municates with its six nearest and eight third-nearest neighbours. The fluid dynamics 
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is, as usual (see Higuera et al. 198£; Ladd 1994) encoded at each site by a distribution 
function fi, where the subscript obeys ^ i ^ 14. This ascribes weights to each of 15 
velocities c^: one null, six of magnitude 1, eight of magnitude with directions such 
that each velocity vector points toward a linked site. In order to model the binary fluid, 



a second set of distribution functions, gi, is also used, following Swift et al. (1996). The 
f's are defined such that 

i 

where the sum is over all directions, i, at a single lattice point, while for <?j the same sum 
gives the order parameter, 

5>i = £ (6.2) 



(At this point, our algorithm departs from that of Swift et al. (1996), who would have 
4>p on the right. The two methods differ only by terms that vanish in the incompressible 
limit of interest, where p — > 1 everywhere.) 

The momentum, pv a (with a a cartesian index) is then given by 



pv a 



where Ci a = (cj) Q . The full pressure tensor, V a /3, is given by 

V a = 22 fiCiaCi/3- 



(6.3) 



(6.4) 



This expression includes not only the conservative stress V^n but also dissipative (viscous) 
contributions, and a trivial 'kinetic pressure' pv a vp which arises in any fluid moving at 
constant velocity v. 

The distribution functions fi,gi obey discrete evolution equations involving simple 
first-order relaxation kinetics toward a pair of equilibrium distributions: 



Mr + Ci ,t + i)- Mr, t) = -{fi - ft 9) )/n 

gi(r + Ci,t+l)-gi(r,t) = ~{g t - g i f q) )lT 2 



(6.5) 



91 ")/T 2 (6.6) 
thus defining two relaxation parameters n , t%. In our use of the code we select t% = 1 

Jeq) 



which causes gi to be reset to g\ each time step. The viscosity is determined by n, with 

r\ = (2ti — l)/6 in lattice units. The equilibrium distributions, /f e<3 ' and g\ eq \ can be 
derived from (6.1) - ( |6.3[ ), along with the condition that the order parameter is advected 



by the fluid, 



i 

and that the pressure tensor and chemical potential at equilibrium obey 

i 

^2 9^ C ia Ci^ = MpS af 3 + (f>V a V0. 



(6.7) 

(6.8) 
(6.9) 



The parameter M controls the order parameter mobility M via MAtfa — 1/2) = M , so 
that M = 2M in our case. (Note that in |Kcndon et al. (1999Q and |Cates et al. (1999[ ), the 
quoted values of M are in fact M values and should therefore be halved to give the true 
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order parameter mobility.) The second term on the right in (3.£) is the trivial 'kinetic 
pressure', with an analagous term in Eq.(6.9). 

By expanding f^ , g^' to second order in velocities and solving for the coefficients 
one obtains 



/| 9J = jA„ + 3u Q C iQ + -V a VpCi a Cii3 - -V + GapCiaCip > . (6.10) 

Here, v is an index that denotes the speed, 0, 1, or -\/3, and lo v , A„ and G Q ^ are constants 
given by 

cj = 2/9; wi = 1/9; w 3 - 1/72, (6.11) 

A = ^ - Jlrp t/l ; A 1 =A 3 = iTrP* h , (6.12) 
2 2 p 



G Q/3 = - ^L T rV th . (6.13) 



9 t-,4/1 3(5 a; 3 TvntA 

The equilibrium distribution for the order parameter, g| ei? ^ , is the same as for /| 



with T 7 replaced by M/i 1 in the above equations. The above results follow Swift et 



(1996 ), generalised to three dimensions. 



To complete the model specification, one must introduce expressions for the pressure 
tensor an d ch emical potential derived from the free energy functional. In this study we 
chose in (2.3) a simple '(/> 4 ' model for V(</>): 



(6.14) 



Ffap] = Jdr j^ 2 + ^ + ^(V^) 2 + iplnpj 

where A < 0. (The term in p is discussed below.) With this choice one finds 4>* = 
±(-A/B) 1 / 2 , and using (fT|), Q, 

a = (-8kA 3 /9B 2 ) 1 / 2 (6.15) 
/i = Acj) + Bcj) 3 - nV 2 (j). (6.16) 

The equilibrium interfacial profile in given by 4>/4>* = tanh(g/£o), where g is the normal 
coordinate introduced previously, and 

Co = {-K/2A) 1 ' 2 (6.17) 

is a measure of the interfacial width. 

An important addition to ( |6.14 ) is the term dependent on density p, here chosen as 
an 'ideal gas' type contribution (up to the factor 1/3). This gives a diagonal term in the 
thermodynamic pressure tensor, which becomes 



v% = jj + ^4> 2 + - K^ 2 cj> - f (V0) 2 | s Q0 + K (d a <t>)(d 0( f>), 



(6.18) 



so that the thermodynamic stress obeys = (p/3)S a p + V^jf m . Thus in practice 
p/3 = P, which is the isotropic pressure contribution normally viewed as a Lagrange 
multiplier for incompressibility. But in fact, our LB algorithm does not know that the 
fluids are meant to be incompressible; instead the ideal-gas term is relied upon to en- 
force incompressibility to within acceptable numerical tolerances. (This avoids a sepa- 
rate calculation, at each time step, of the fluid pressure P and renders the algorithm 
local.) Compressibility errors can be minimized by increasing the coefficient of plnp in 
( |6.14 ), which would require a shortened time step, or for fixed time step, by reducing 



14 V. M. Kendon, M. E. Cotes, J-C. Desplat, I. Pagonabarraga and P. Bladon 



the magnitudes of A, B, and n together so that an acceptable level of incompressibility 
is maintained. The second route is followed here. 



The chosen functional, (S.14), has a number of points in its favour for numerical sim- 
ulation. The main terms in V^a and p are simple powers of 0, so are easy and quick to 
evaluate. (Models involving logarithms or trigonometric functions (Swift et al. 1996) pay 
a heavy price in computational efficiency.) Further, the shape of the "0 4 " potential is 
fairly smooth, avoiding very steep gradients that might lead to inaccuracy and instabil- 
ity when approximated numerically on a lattice. We nonetheless need to evaluate spatial 
gradients of </>; this is done using all 26 (first, second and third) nearest neighbours on 
the D3Q15 lattice, so that numerically 



<9 Q </>(x) 



(r + Ci) 




) - 26^(r) 



(6.19) 



(6.20) 



Note that these are not the only possible choices and that others will only coincide when 
all gradients are small on the scale of the lattice. There may be considerable scope for 
improvement by optimising the choices made, but we leave this for future work. 

As usual in LB, we have chosen a lattice with enough symmetry to ensure that the 
rotational invariance of the fluid mechanics is faithfully represented (Higuera et al. 1989). 
However, this does not guarantee that the same holds for the thermodyamic properties. 
With our choice of free energy functional, (6.14) and the above gradient discretisation, 
rotational invariance in the thermodynamic sector is recovered only when all order pa- 
rameter gradients are weak, which in principle requires £o ^ 1- In practice, a compromise 
is necessary; we return to this in § |3.2| below. 

Finally, the hydrodynamic behaviour of the LB technique requires detailed comment. 
The hydrodynamic equations that correspond to LB can be obtained by making a 
Chapman-Enskog expansion of the Boltzmann equations (3.5 6.6). If we consider the 



expansion for the distribution function / (which relates to the fluid momentum) we 
arrive at 



dr(pv a ) + d a (pv a vp) = 



d 



V ( dpv a + d a v/3 - -SapdyVy ) + S,d 7 v 7 S al 3 



[v a d 7 V^ em + v^v c J; em 



3 a/ 



f3 /3 9 7 (pw a w / 3U 7 )(6.21) 



p • - — * P 

where rj and £ are the shear and second viscosities, respectively. For our single-relaxation 
LB scheme £ = 2^/3. 



The first line of Equation 3.21 corresponds to the standard Navier Stokes equation, 
and shows that, through these terms, the model recovers both the compressible and 
incompressible features of isothermal hydrodynamics. (J 

The second line in Equation 3.21 contains spurious terms, which arise partly because 
the enthalpic interactions that lead to the non-ideal behaviour of the LB fluid (or fluid 
mixture) are introduced through equilibrium information only. (In a Hamiltonian system, 
the same interactions that perturb the equilibrium state away from an ideal gas would 
also be responsible for the dynamics.) Of the two terms, the second one is not Galilean 

f Note that when modeling a macroscopic length, L, the ratio ri/(c 3 pL) will be larger than 
in typical real fluids; due to the presence of the lattice, LB does not have a sharp separation 
between the c ompressible and inc ompressible time scales. In this respect, it resembles a high 
viscosity fluid, Hagen et al. (1997). 
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invariant but is cubic in the velocity. It will be negligible for small velocities (recall 
that the LB algorithm anyway requires fluid velocities that are small in lattice units). 
The first term can be decomposed into a Galilean-invariant part, which is a product 
of gradients of the pressure and velocity, and a non-Galilean invariant contribution. 
The product term is small compared to the Navier-Stokes terms, which are linear in 
such gradients, whenever these are weak; the non-Galilean invariant term is small under 
the same conditions. There have been recent proposals to enforce Galilean invariance 



within compressible multiphase lattice-Boltzmann schemes, Holdych et al. (1998), and 
this is a desirable feature of future algorithms. Nonetheless, under circumstances where 
all hydrodynamic field s var y smoothly on the lattice scale, the spurious terms appearing 
in the second line of ( 3.21 ) will typically be smaller than the retained ones of the NSE 
equation in the first line. 

The evolution equation for the order parameter can be obtained analogously, by per- 



forming a gradient expansion of the linearized Boltzmann equation (6.6). This leads to 



4> + d a ((t)v a ) 



T'2 



yth 

13a 



(6.22) 



Equation ( |6.22 ) has the usual form of a convection-diffusion equation, so long as one 
chooses M as a constant, except for the last term. As with (6.21), this Galilean-invariant 



term arises as a result of the way in which the non-ideality of the fluid mixture is in- 
troduced; like the first spurious term in (6.21) it contains one higher derivative than 



the term d a V^p that enters the Navier Stokes equation, and is expected to be small for 
similar reasons. In Appendix ^ we confirm explicitly that, in the incompressible limit 
only, this spurious term does not modify the hydrodynamic modes of a binary mixture, 
at the level of a linearised expansion about a uniform quiescent fluid. 

In summary, for a nominally incompressible fluid the correct fluid behaviour is re- 
covered in the only regime where it can justifiably be expected, namely, when all the 
hydrodynamic fields vary slowly on the lattice length scale. Within the current LB al- 
gorithm, one also depends on having only slight fluid compressibility: this eliminates a 
spurious coupling between order par ameter and momentum fluxes (Appendix |A|, and 
Swift et al. 1996 ; Holdych et al. 199S| ). Because of this it is important, with the current 
algorithm, to monitor closely the actual fluid flow. 



7. Parameter Steering 

Since it is possible in a single simulation to sample only a small piece of the l(t) 
curve, it is necessary to work one's way along this curve via a series of different runs. 
This means varying Lq relative to the lattice constant so that an appropriate window of 
reduced length scales I lies within the range of the simulation. Put differently, one must 
set, in lattice units, 1 <C Lol(t) <C A. 

In the simulations, however, we need to choose not one parameter value (Lq) but seven 
(A, B, k, p, M, and also tj., tq). Few of these parameters can, in practice, be set indepen- 
dently: an unguided choice would typically produce a simulation that either did nothing 
over sensible timescales or became unstable very rapidly. To avoid these outcomes, care- 
ful parameter steering is required. This was done by semi-empirical testing using small 
(workstation) simulation runs until satisfactory choices emerg ed. Runs on 96 3 and 128 3 
lattices were used to confirm these before committing the resources for 256 3 runs. The re- 
sulting parameters are summarized in tables ^|and|^. The guiding principles that emerged 
from this process, as well as the actual parameter values, are of some interest to those 
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planning future work with LB. They are now summarised briefly, followed by a discussion 
(§ ||) of several validation exercises that were then undertaken. 

In essence, our navigation of the l(t) curve involves steering three parameters (77, ex, M) 
at fixed values of the remaining four [B/A, k/A, /?, T2). Firstly, the (mean) density p can be 
set to unity without loss of generality; we do this. Second, we set <f>* — 1 by choosing B = 
—A; in terms of the original definition of the order parameter (0 = (ha — tlb) / '{tia + tib)) 
this amounts to a simple rescaling of <f> — ► (f>/4>*. Varying k in proportion to A then gives 



3.2) 



control over the interfacial tension a (6.15) while retaining a fixed interfacial width £0 in 



lattice units (Eq. (6.17)); this keeps thermodynamic lattice anisotropics under control (i 



To achieve an efficient simulation, one requires the interfacial velocity L to be of order 
0.01 in lattice units during the main part of each run. (Any slower will exhaust resources; 
any faster will give compressible and inaccurate fluid motion, and, in all likelihood, 
numerical instability.) At each point on the l(t) curve, this gives, a posteriori a relation 
between ij and a. Thus to access large I one clearly requires small L — rj 2 /(pa); but 
to avoid compressibility problems (§ B^J) this must be done by reducing viscosity rather 



than increasing interfacial tension. Maintenance of numerical stability (§ 8.1) requires in 
fact that we decrease a with decreasing r\\ however, these factors do not cancel in Lq and 
a wide range of values (about six decades) can stably be achieved. Thus we were able to 
explore the viscous, crossover, and inertial regimes; these various regimes are delineated 



quantitatively in § 9.2 below. 

Setting the correct mobility M is crucial throughout. Across the whole l(t) curve, one 
has to ensure that M is large enough that interfaces relax to local equilibrium on a time 
scale fast compared to their translational motion. But if M is made too large, residual 
diffusion becomes a significant contributor to the coarsening rate, contaminating the 
data. This tradeoff can be eased in principle by going to larger system sizes than those 
currently available. It could also be improved by making M a function of 4>, setting 
(for example) M = M (l — <fi 2 ). This would have the effect of giving strong diffusion 
only where it is needed, in the interfacial region. However, implementation of this within 



LB is not trivial (Swift et al. 1996); specifically it is not enough to make M in (6.9) 
0-dependent|. 

It is not surprising that mobility is a limiting factor at large L (viscous regime, small 
I): diffusion will always enter if the fluid flow is slow enough (high enough rj). But mobility 
factors also come into play at the inertial end (small L , large I): in physical units, the 
interface in this regime is unnaturally wide and to maintain it in diffusive equilibrium 
(and keep the algorithm stable) again requires relatively large M. These cause residual 
diffusion which, for our system sizes, limits from above the range of I accessible. 

Finally, we found that accuracy in the viscous regime (small I, large L ) is compromised 
when the viscosity becomes too large (of order unity, in lattice units). The signature of 



this is an apparent breakdown of energy conservation (see § 11. 5| ) . We are not sure of its 



origins, but note that too large a viscosity causes the dynamics of momentum diffusion 
and that of sound propagation (density equilibration) to mix locally. In an almost steady 
flows this should not matter, but in the small I regime the viscous and interfacial terms 
in NSE are both numerically large. In principle these balance to give negligible fluid 
acceleration but their numerical cancellation may be imperfect. Although any such local 
accelerations are numerical in origin, the response to them may need to be accurate, if 



f Since in ( EL22 ) M enters in the form V 2 (M/^) rather than as V.(MV(i) as would be 
required if M were not constant. 
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0.0021 
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Table 2. Parameters used in 256 3 lattice-Boltzmann runs. 
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Table 3. Parameters used in 128 3 lattice-Boltzmann runs. 



the global physics is to be handled correctly. We speculate that this is a limiting factor 
in our exploration of the lit) curve at the lower end. 

In this study, the largest system size was A 3 = 256 3 , although due to disk storage lim- 
itations, the results from this system size were analysed only after coarse-graining down 
to 128 3 . The coarse-graining was done by averaging over blocks of eight neighbouring 
lattice sites to create one coarse-grained value. Runs at 128 3 and 96 3 were also done, and 
results for all calculated quantities were compared between 256 3 and 128 3 runs with the 
same parameters, to identify any effects of coarse-graining. The main 96 3 and 128 3 , 256 3 
simulations were run respectively on the EPCC Hitachi SR-2201 machine (4 processors) 
and the EPCC Cray T3D (64 and 256 processors). Follow-up studies used in some of the 
velocity analysis work, and for additional visualisation, were made on the CSAR Cray 
T3E at Manchester. 

Typical runs required, in the 256 3 case, around 3000 T3D processor hours CPU, and 10 4 



time steps to reach the point where finite size effects set in (see § 8.5). All simulation were 
run with periodic boundary conditions; the initial configuration was always a completely 
mixed state, with small random fluctuations. For each run, the order parameter, 4>, and 
the fluid velocity vector at each lattice site were saved periodically for later analysis. The 
sampling frequency was limited by the available disk space. Typically, data was saved 
every 300 time steps giving, over a run of 10 time steps, around 4Gb of data. 



8. Validation and Error Estimates 
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Figure 3. Interface profiles (<j>) and gradient profile (V(f>), for spherical droplet equilibrated in 
the opposite fluid. Data has been collected in bins of width 0.1 lattice spacings, and the error 
bars are one standard deviation. The sphere radius is 32. Left: interface width set by £o = 0.57. 
The theoretical profile is tanh(g/£n), for a flat interface, where g is a coordinate normal to the 
interface. Right: interface width set by £o = 0.88. 



8.1. Numerical stability 

The LB method is not generally stable. In fact, our experience suggests that, whatever 
parameters are chosen, any run would eventually become unstable if continued for long 



enoug h; this is not dissimilar to some molecular dynamics algorithms, Allen fc Tildcslcy 



(1987 ). During testing, a reliable picture was acquired of the characteristic way in which 
this happens. When the inaccuracies have built up to the point of failure, the velocities 
become very large over a small number of time steps until numerical overflow causes the 
code to stop running. There seems to be no danger of taking data from a period when 
the system might be far from accurate but still apparently running successfully, since the 
onset is so rapid. Thus there are several runs among the set used for final data analysis 
where the run ended prematurely due to instabilities, but the data prior to the instability 
has been considered sufficiently reliable to be used. 

8.2. Anisotropy and interfacial tension 

The elimination of lattice anisotropy in the thermodynamic sector of the model requires 
£o 3> 1 in lattice units, to ensure that the interfacial tension a is independent of interface 
orientation. In practice this goal must be balanced against other demands. To test the 
extent of the problem, a spherical droplet (radius 32 lattice units) of fluid B surrounded by 
fluid A was allowed to equilibrate. The interface profile was then measured by evaluating 
the mean and standard deviation of the order parameter <f)(r) at various radii r (binned 
on the scale of 0.1 lattice units) from the droplet centre. The result is presented in figure 
H for £ = 0.57 and for £ = 0.88. Note that the 'width' of the interface, as judged by 
eye, is actually about 5£o- 

A closer look at the droplet shape (not shown) in each case reveals that the sphere 
has deformed slightly by squeezing along the Cartesian lattice directions and expanding 
along the diagonals. This deformation is about 3.5% for the narrower profile and about 
1.5% for the wider profile. A similar test, done with a sphere of radius 31.5, confirmed 
that this deformation was not due to any tendency of the interface to lock onto specific 
lattice sites but purely from anisotropy of the tension. 

If all else were equal, the wider interface would be chosen. However, the computational 
penalty for wider interfaces is severe. To maintain these in local equilibrium, the mobility 
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—A, B k n Ma theory a measured 

0.083 0.053 1.41 0.1 0.063 0.055 

0.063 0.04 0.5 0.5 0.047 0.042 

0.0063 0.004 0.025 4.0 0.0047 0.0042 

0.0031 0.002 0.0014 8.0 0.0024 0.0021 

0.0013 0.0008 0.0005 10.0 0.00094 0.00083 

Table 4. Interfacial tension, theoretical and measured values. 



M must be high enough to allow diffusion across several £o on a timescale faster than 
fluid motion. For the wider interface (£o = 0.88) the resulting residual diffusion then 
contaminates much of the remaining L range. It was thus found necessary to sacrifice 
some isotropy for efficiency, and the narrower profile with £ = 0.57 was used for the 
main runs in this work. The resulting anisotropies are marginally detectable by eye in 
visualisations of the interface for the spinodal system (e.g. figure [To| below). We estimate 
that they contribute systematic errors of a few percent to the growth rate L{T) 1 which 
is comparable to other sources of error. 

The mean interfacial tension was measured for each parameter set by allowing an 
interface to come to equilibrium and numerically performing the integration in ( |2.7| ). 
Both terms were evaluated, and an average taken over various configurations. This gives 
values for the interfacial tension, shown in table [|, that are systematically about 10-15% 
smaller than the theoretical values. (The statistical error is a few percent.) The difference 
is due to the narrow interface leading to inaccuracies in the gradient calculations. But 
as far as the simulation is concerned, this systematic effect is removed by our using the 
measured value of the interfacial tension in subsequent calculations of L$ and To. 

8.3. Local equilibrium and residual diffusion 

Errors in the intereface-driven dynamics can arise if the interface is not maintained 
in local equilibrium. This was tested as follows. Since the bulk fluid is fully separated 
((f) = ±1), one expects 1 — (\<f>\) oc A/V oc 1/L where A/V is the area per unit volume and 
angle brackets are a real-space (site) average. Within a given run, any departure from 
constancy of the product L(l — (\4>\)) is thus an indicator that the interfaces are failing to 
keep up with the evolution of the surrounding fluid. (This product could have different 
asymptotic values in the viscous and inertial regimes, so the product need not be the 
same in different runs.) At the lowest values used for the mobility M (deep in the viscous 
regime) there was measurable deviation from constancy, from which the nonequilibrium 
deviations in a were estimated to be of order 5%. Any deviations in the inertial regime 
were, however, smaller than this. 

Careful checks were made to exclude residual diffusive contributions to the coarsening 
process. This was done using comparator runs in which the viscosity was set to an 
extremely large value so that coarsening was purely diffusive. (Such runs are depicted 
in figure ^ below.) From this, the diffusive coarsening rate was found as a function of 
domain size. Then for the full run (with fluid motion reinstated) all data was excluded for 
which this diffusive coarsening rate exceeded 2% of the full rate. This whole procedure 
was repeated with a limit of 1% instead of 2% on the residual diffusion. The values of 
the fitted exponent a as per ( |4.9| ) (given in the last column in table ||), did not change 
beyond the estimated errors so the limit of 2% diffusion was taken to provide sufficient 
accuracy. 
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Figure 4. Left: ratio of the radial to transverse velocity components in Fourier space for various 
runs. Right: velocity structure factor showing relative magnitude of the Fourier space velocity 
at different wave vectors. Wavevector axis is labelled by wavelength in lattice units. 



The result of this choice was exclusion of data with L < L m i n ~ 15 — 25 (varying 



somewhat between runs). Had a wider interface been used (see § 3.2) then by the same 
criterion L m i n would be much larger giving very little usable data. 

8.4. Compressibility and small scale structure 

The Ludwig code will only work correctly at low Mach number. This requires L <C c 
where the sound speed is c = in lattice units. Since in our simulations L is of order 

0.01, we expect our the binary fluid mixture to remain incompressible (V.v = 0), at least 
at length scales larger than a few lattice sites; in Fourier space, we expect k.v(k) = 
at all but high k. figure ^ shows the RMS ratio of the radial to the transverse veloc- 
ity components in Fourier space as a function of wavenumber, and also the spherically 
averaged velocity structure factor, S v (k) = (|v(k)| 2 ), for various runs. Also shown for 
comparison is single fluid turbulencejj], g enerated using pseudo-spectral direct numerical 
simulation (DNS) code by Young (1999| ), and a LB run with a single fluid (no interface) 



but otherwise the same parameters as Run031 (inertial region). 

At low wavenumbers the sytem is incompressible. At higher wavenumbers, there is 
some compressibility, whose effect varies in the different growth regimes. In the viscous 
regime, the longitudinal/transverse ratio rises with fc, but the velocity structure factor 
shows that that all velocity components become small at high k and contribute little to 
the overall dynamics. This is still true in the crossover region, where the compressibility 
ratio is highest; a peak in S v (k) is found at a wavelengths around 3 lattice spacings. In 
the inertial region, this peak shrinks, and splits into two (at around 3.5 and 2.5 lattice 
spacings). The transverse velocity component is now larger although still an order of 
magnitude smaller than the velocity at the peak of S v (k). 

Comparison with the single fluid turbulence, as simulated by both DNS and LB, shows 
that these peaks in S v (k) are mainly due to the presence of the interface. Their presence 
only in the crossover and inertial runs suggests that perhaps capillary waves are forming 
on the interface giving structure the velocity field on scales of the order of the interface 
width. Subsequent visualisation work showed that underdamped wavelike motion of in- 
terface is undoubtedly present at large I, Desplat et al. (2000fi| ), but predominantly at 



f The single fluid turbulence simulation method sets the radial component identically to zero 
thus guaranteeing perfect incompressibility. 
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wavelengths much larger than the interfacial width. Another argument against the capil- 
lary wave explanation is that no similar bumps are seen in the order parameter structure 
factor S{k) (figure |). 

The nature of the velo city fields close to the interface cer tainly deserves further investi- 
gation (see, for example, Theissen, Gomppcr fc Kroll 1998 , for related work on a different 
system). Meanwhile, to have some compressibility on the length scale of the interface it- 
self appears unavoidable within current LB. Specifically, in the immediate vicinity of the 
interface the various diagonal terms in the chemical contribution to the pressure tensor, 
(2.10), are individually large, although these should nearly cancel for a slowly moving, 
weakly curved interface. Any numerical error here will lead to local deviations in the 
fluid density p, even if the bulk fluid motion is effectively incompressible everywhere else. 
On molecular physics grounds also, some coupling between density and order parameter 
can be expected at the interface between otherwise incompressible fluids. Such coupling 
is present in real physical systems, but care is needed with the current LB code where 
compressibility effects also bring violations of Galilean invariance (§ pf). 



8.5. Finite size effects 

Various estimates were made of when our (periodic) boundary conditions started to sig- 
nificantly influence the behaviour of L(T). This included several comparisons of different 
sized runs with the same values for other simulation parameters. On this basis, the data 
for the 96 3 and 128 3 was pruned at L = L max = A/4 before analysis, and the 256 3 
runs terminated at this point. This criterion is much more conservative than in some 
previous work (e.g. Jury et aL 1999& ), and, given L m i n , limits the range of L accessible 
in a single large run to about half a decade. To balance this, averaging over different runs 
with the same parameter values should not then be necessary, since one has in effect 
(A/L max ) 3 — 64 different (albeit correlated) samples being simulated within each run. 
Indeed, in the crossover and the inertial regime, we saw no sign of statistical fluctuations 
in the L(T) plots. 

Interestingly, the same was not true for the extremely viscous runs, which showed 
somewhat erratic statistics (see § ^|). One possible reason for this is the presence of 
correlations, in the velocity field, over much larger length scales than L(T), causing the 
local coarsening rates in different parts of the simulation to fluctuate coherently. Long 
range velocity correlations are, in fact, clearly visible in the structure factor S v (k) shown 
in figure ||. Specifically, for the most viscous run analysed (Run 027, L = 36), S v (k) 
shows no sign of saturating at low k; instead the data suggests a power law divergence, 
and is consistent with S v (k) ~ k~ 2 . (A theoretical argument leading to this result for 



the viscous regime is given in § 10.1.) In real space this translates into a long range, 



\jr velocity correlation extending to either the system size (which is the likely case in 
any simulation) or some large physical length scale beyond which the purely viscous 
approximation (Stokes flow) breaks down. 

If this is correct, it could be practically impossible to avoid finite size effects when sim- 
ulating the viscous regime. The most benign outcome is if the main effect is to correlate 
(rather than alter) local coarsening rates; this could be countered by averaging over a 
number of different runs (Jury et al. 19996; Laradji et al. 1996). However, this would 
have to be done for several system sizes before concluding that no other finite size effects 
were present. 
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Figure 5. Left: structure factor, S(k), for Run028 (viscous regime) for timesteps 14000 - 19000, 
L(T) = 38 - 52. Right: S(k) for Run032 (inertial regime) for timesteps 11000 - 17000, L(T) = 45 
- 64. Different open symbols denote different times T; filled circles show the same data corrected 
at low k for discretisation effects (see text). 



9. Order parameter results 

We now present our results for the time evolution of the interfacial structure. These 
results can be extracted directly from knowledge of the order parameter using well- 



established procedures ( see Jury et al. 19996 ; Appert et al. 1995 ; Laradji et al. 1996 ; 
Bastca fc Lcbowitz 1997| ). We defer to § |10| our explicit analysis of the fluid velocity field. 



9.1. Structure factor scaling 

The first step in the analysis of the order parameter data was calculation of the structure 
factor. The (f> field saved from the simulation runs was processed through numerical 
Fourier transform routines, and the structure factor calculated as: 



S(k) = — V 0(k)0(-k), (9.1) 

fc~7r/A<|k|</c+7r/A 

where 0(k) is the Fourier transform of the order parameter, and is the (actual) number 



of lattice sites in a shell of radius k and thickness 27r/A in Fourier space (compare (4.2)). 

Dynamical scaling requires that, in reduced physical units, not only the characteristic 
length l(t) but also the statistical distribution of different interfacial structures should be 
the same for each I. In either the viscous or the inertial regime, therefore, the structure 
factor S(k) should asymptotically collapse onto a single plot when appropriately scaled, 
so that in simulation units 

L- 3 S(k) = f(kL), (9.2) 

with a different function f(kL) in each of the two limits. (More generally, dynamical scal- 
ing allows L~ 3 S(k) — f(kL,l), so that the viscous and inertial asymptotes are f(kL,0) 
and f(kL, oo) respectively.) figure || shows plots of S(k) scaled in this way for Run028 
and Run032, representative of the viscous and inertial regimes respectively. 

The collapse of the structure factor data within each run is good (figure ||) for length 
scales larger than about twice the interface width. (The latter is marked as 2£ on the 
graphs, with £ = 5£o-) With our definition of L, the peak occurs at kL just less than one. 
To the right of the peak there is a shoulder, followed by a reasonable approximation to a 
k~ A Porod tail. (The Porod tail represents scattering from a weakly curved interface and 
should be found in the region i^k- 1 < L(T), see |Bray (1994J ); but between £ ~ 3 and 
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Lmax — 64 there is barely room to observe it cleanly.) The ragged sections of S(k) in the 
low kL region corresponds to the first two fc-shells which have few k points and so poor 
statistics. The filled symbols are the same data corrected to allow for the fact that the 
average value of |k| in such a shell differs from the nominal shell radius; the corrected 
result suggests no deviation from scaling even at low k, although the data there is less 
reliable than in the high wavenumber region. 

The collapse between different runs (not shown) is also good, so long as one compares 
runs chosen within either the viscous or the inertial regime. However, as is visible from 
figure H, the shape of S(k) does evolve significantly between one regime and the other. 
In particular, the shoulder to the right of the peak is lower in the viscous regime than 
the inertial regime. This implies that the domains are a subtly different shape in real 
space, perhaps more evenly rounded in the linear regime since the peak is effectively a 
little sharper. This may be linked to an increased number of relatively narrow necks in 



the inertial runs (large I), as first suggested by Jury (1999 ) and recently confirmed by 
direct visualisation of LB data, Desplat et al. (2000 1|). Our st r ucture factor results, taken 



piecewise, are compatible with those of Jury et al. (1999a ), Appert et al. (1995 ), and 



several other authors (see Appendix |BJ). However, our study is the first to cover a wide 
enough parameter range to show a clear distinction, in the shape of S(k), between the 
viscous and inertial regime. 

Runs in the crossover region also show reasonable data collapse within each run, with 
a shape intermediate between the two shown in figure ^, and very similar to that found 



by Jury et al. (1999a) in the same region of the l(t) curve. Note that a good collapse, 
within or between runs, cannot be expected a priori in the crossover region. It arises 
because the /-dependent scaling function f(kL, I) in fact evolves so slowly with I that 
any data spanning less than a decade or two in I is insensitive to the I dependence. This 
is a consequence of the extreme breadth of the crossover region (quantified below) . 

9.2. Evolution of the characteristic length scale 
The characteristic length scale L(T), defined via ( |4.l| ), has been calculated for the eight 
256 3 runs in table |[ The order parameter data was coarse-grained to 128 3 before analysis, 
but comparison with smaller runs confirmed that there was no effect of this on L(T) 
within the 'good data' range. The latter is defined as L m i n < L < L max , with L m i n fixed 



by our criterion on residual diffusion (§ 8.3) and L max = 64 as required to exclude finite 
size effects (§ |3~5| ). Figure ^| illustrates how the fitting was done. 

To parameterize the time dependence of L(T), the 'good data' was fitted, for each run 
separately, to the following form 

L = v(T-T int ) a , (9.3) 

(equivalent to ( f4.9| )) where v, Ti nt and a are fitting parameters. A nonlinear curve- 
fitting utility was used to create the fits, which all fell within a specified tolerance of 
1%. However, some trade-off is possible between the three fit parameters and a realistic 
uncertainty estimate for the exponent, a, is around 10% for the first three runs in table 
||, and 5% for the rest. The fits are shown in figures ^ and 0, which also shows the 



diffusion-only data used to determine L m i n as described in § 8.3. The fitted results are 
summarised in table [|. 

The data show a values ranging from 1.12 to 0.66 with a decreasing trend as Lq is 
decreased. Certainly, an increasingly negative curvature of the L(T) plots with decreasing 
Lq is apparent from figures [?], |^. However, the resulting fit parameters were relatively 
erratic for the three runs of largest Lq (expected to lie in the viscous regime). Indeed, we 
found a = 0.86, v = 0.023 and a — 1.16, v — 0.0012 for two runs with the same nominal 
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Figure 6. LvsT graph (unsealed) for Run030, illustrating the fitting procedures. 
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, for 256 3 runs 


The parameter vd is the fit parameter 



corresponding to v in the presence of diffusive growth only (see § p.3| 



Lq. This was partly due to a relatively ill-conditioned fit as can be appreciated from 
figure |t]. (A secon d po ssible cause of the erratic fits is the presence of long-range velocity 
fluctuations; see § |3.5[ ) 

Therefore it was decided to refit the data for the three most viscous runs, imposing 
a = 1, the anticipated value. This yielded much better consistency among the fitted values 
of v, which with viscous scaling should obey (To/Lq)v = bx, where b\ is universal; with 
the forced linear fits this was indeed the case with b\ extracted as 0.073, 0.072, 0.072±0.02 
for the three runs under discussion. Subject to this, we obtain a range of values of a from 
1.0 (Run028, Run022 and Run033) to 0.67 (Run019), with intermediate exponents 0.95 
(Run029), 0.80 (Run020) and 0.75 (Run030) at intermediate L . This suggests that the 
simulations have indeed covered the viscous, crossover and inertial regions. However, the 
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Figure 7. Fitting L(T) and L D (T). Upper left: Run028, L = 36. Upper right: Run022, 
Lo = 5.9. Lower left: Run033, Lo = 5.9. Lower right: Run029, Lo = 0.95. Solid lines indi- 
cate full set of recorded L(T) data, + indicates data points used for fits with L min set by 2% 
diffusion, Q indicates data points used for fits with L m i„ set by 1% diffusion, A indicates data 
points used for fits to diffusion- only data, table o summarises the main fit results. 



ultimate test of this is to convert to reduced physical units and construct the l(t) curve 
explicitly. 

9.3. Universal scaling plot for l(t) 
Our method for combining the data from different simulation runs to give the l(t) curve 



follows Jury et al. (19996). As apparent from definitions (1.1), the only fit parameter 
that is actually needed when converting L(T) data to reduced physical units {1(f)) is 
the intercept, Ti„ t . Then one uses the known density and viscosity, and the measured 
interfacial tension, to complete the conversion. 

Figure || shows the 1(f) data from all the runs in table || combined on a single log-log 
plot. Note that for the two runs of L = 5.9 the resulting data collapse is much improved 
by the forced linear fit (giving two very different values of the nonuniversal offset T int 
instead of two disparate values of a). With the former, the two datasets overlie on the 
1(f) plot but with the latter they do not; this helps to vindicate our choice of fit. Apart 
from a similar reservation about force fitting a for the most viscous run (Lo = 36), the 
1(f) curve is free of adjustable parameters. Although we did not have resources to cover 
the entire curve with data, there is no evidence for any breakdown of universality: the 
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Figure 8. Fitting L(T) and L D (T). Upper left:Run020, L = 0.15. Upper right: Run030, 
L = 0.01. Lower left:Run019, L = 0.00095. Lower right: Run032, L = 0.0003. Solid lines 
indicate full set of recorded L(T) data, + indicates data points used for fits with L min set by 2% 
diffusion, Q indicates data points used for fits with L m i„ set by 1% diffusion, A indicates data 
points used for fits to diffusion- only data, table o summarises the main fit results. 



various runs do appear to lie on a smooth underlying curve. (In particular, the two most 
inertial runs virtually join up.) 

The apparently universal l(t) curve shows scaling that is first linear (I = bit, with bi = 
0.072 ± 0.02), then passes through a broad crossover region before reaching I = b 2 /^t 2 / 3 
(with 62/3 = 1.0 ± 0.05) at large I, t. The positions of the crossover and inertial runs 
on the graph are in keeping with the trend for the scaling exponent, a, fitted directly 
from each run. This confirms that the exponents determined from our fitting procedure 
do accurately reflect what is going on in these simulations. The extreme breadth of the 
crossover regime, 10 2 < t < 10 6 , justifies the use of a single exponent to fit each run (4.9), 
( |9.3| ) even in the crossover: no single run is long enough to see a change in exponent from 
beginning to end beyond the estimated errors. Th ere is no hint that the exponent is 
reducing still further to a ^ 1/2, as predicted by Grant fc Elder (1999 ), although a 
further crossover beyond the range of l,t reached in these simulations cannot be ruled 
out. 

Recall that intersection of asymptotes on the l(t) plot defines t* , the characteristic 
crossover time from viscous to inertial behaviour. As mentioned previously (§ 0), scaling 
theory says only that t* is 'of order unity'. The measured value is close to lCr, a value 
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Figure 9. Scaling plot in reduced variables (L/Lo, T /To) for 256 3 LB data. Dots (left to right) 
are the runs in table ^ (top to bottom). Dashed lines show free exponent fits for the first three 
data sets in table ^ for comparison with linear fits (dots) . 



that should raise no eyebrows in the turbulence community but may do so among workers 
in phase separation kinetics. Since 6 2 /3 is very close to 1, the largeness of t* can be traced 
to the smallness of b± and to the relatively minor change in exponent on crossing from 
viscous to inertial scaling: for by its definition, t* = {bi/b v ) 1 ^ av ^ a ^ where subscripts i, v 
signify viscous and inertial values. 

Note too, the huge range of scales covered by the combination of eight simulation runs: 
five decades of length and seven decades of time. This achievement is only possible by 
fully exploiting the expected scaling. This means that, although our work is capable of 
falsifying the scaling hypothesis (our l(t) plot might not have joined up, and might yet 
not do so when more data is added), its non-falsification in our work may not represent 
persuasive proof that the scaling is true. 

For, as mentioned previously (§ ^), to navigate the l{f) curve we are forced to correlate 
the simulation parameters in a systematic way. Hence if the coarsening rate was in fact 
dependent on M, say (for example by being pinchoff-limitcd, Jury et al. (19996)), this 
would not necessarily show up as bad data collapse in figure |9j, since M is strongly 
correlated with Lq and/or To. In principle, however, our parameter steering has no effect 
on data within a run, so that any 'steering-induced' data collapse could in principle 
be detected because curves would not quite line up with their neighbours on the plot, 
although their midpoints would lie on a smooth curve ( [Jury et al. 1999 b\ ). Although we 
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believe this is not happening for our l(t) data, we do not have enough results to entirely 
rule it out, especially as something similar does occur in our own velocity derivative data 
(§0). 

Our results for the evolution of the interfacial structure are compared with those of 
previous authors in Appendix M. 



10. Results for the fluid velocity 

Due to data storage limitations for the largest (256 3 ) runs, our velocity analysis also 
made extensive use the 128 3 runs listed in table ||. The velocity field was analysed as a 
single, continuous field, filling the whole simulation; there is no explicit information about 
the location of the interface between the two fluid phases. However, visualisation of the 
velocity field was done using the AVS package and examples (one viscous, one inertial) are 
shown in figure [l^, where the flow patterns can be compared with the domain structure 
defined by the interface. 

There are almost no prior data on fluid mixtures with which to compare these results. 



A simulation using a pseudospectral method written by Young (f999 ) was therefore used 
to generate a velocity field for single fluid, freely decaying turbulence with a similar 
Reynolds number to those of the spinodal system in the inertial regime. A velocity map 
for this single fluid turbulence is also pictured in figure |h]. 

A trend from locally laminar flow to more chaotic motion is apparent in passing from 
the viscous to the inertial regime. The vorticity map in the latter case is comparable to 
the one for the turbulent single fluid (not shown). However, the comparison is hindered 
by the fact that the interfacial motion at length scale L introduces a 'whorly' velocity 
pattern even in the purely viscous flow regime. A better discriminator between the two 
regimes, pursued elsewhere, comes from watching the time evolution of the interfacial 



structure itself, which is clearly underdamped in the inertial case, Desplat et al. (20006 ). 

fO.f. Velocity structure factor 



The velocity structure factor was introduced in § S.4. For numerical purposes we define 
it (following (fyj)) as 

^W = ^ E v ( k )- v (-k). (10.1) 

The results were shown already in figure [|, where S v {k) is depicted for three of the runs 
in table |[ alongside two calculations (LB and spectral) for single fluid turbulence. These 
structure factors are in unsealed units but in each case correspond to a point during the 
run where the domain size is around 30 lattice units. 



The bumps on the S v (k) curves at high k were discussed in § 3.4. But even apart from 
these, the velocity structure factors have very different shapes in the viscous (Run027), 
crossover (Run018), and inertial (Run031) regimes; these differences are much larger than 
for the order parameter structure factor (figure^). In other words, the geometry of the 
fluid flow is changing much more significantly, as one moves along the l(t) curve, than 
the geometry of the interface. 



We return to this in § 10. 2[ but first address an issue raised in § |8.5| , which is the 



apparent k divergence in S v (k) at low k in the viscous regime. This can be qualitati vely 



explained as follows. In a purely viscous approximation (Stokes flow) the NSE (3.2) 
becomes in Fourier space 

r/fc 2 v(k) = -ik.V th (k). (10.2) 
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Figure 10. Pictures of interface and velocity maps in viscous and inertial regimes. Top row: 
interface (left) and velocity field (right) for Run027 in the viscous region at time step 8500, when 
the domain size is about 21 lattice units. Middle row: interface (left) and velocity field (right) for 
Run031 in the inertial region, for time step 5000 when the domain size is about 22 lattice units. 
Only a 32 3 section of the simulation is shown (the same section for both interface and velocity). 
The velocity is shown for the front two lattice planes only for clarity. Arrow colours indicate 
speed (redder is faster). Bottom row: (left) interface (blue) for Run031 with the red interface at 
a contour of 50% of the maximum vorticity and thus enclosing regions of high vorticity. Velocity 
from single fluid turbulence (right). Turbulence has Re a — 10, matching that of the most inertial 
run, Run031. The vorticity for single fluid turbulence (not shown) is qualitatively similar to that 
shown for Run031 (left). 
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Figure 11. Left: L V (T) (dashed) compared with L{T) (bold) for runs in table || Right: For runs 
in table ^, reduced interface velocity, £(£) (bold) compared with the RMS fluid velocity measured 
in reduced physical units (dotted). 

Here T th contains the chemical term 'p chem which is mainly localised on the interface 
between the two fluids. We now argue that this term is strongly correlated at length 
scales up to the domain size L but not larger ones: this means that, for the purposes 
of low wavenumbers (k -C tt/L) it is a random variable with short range correlation. 
Ignoring for simplicity all tensor indices, one thus has (IT 7 (k)| 2 ) — > \i a constant, as 
fc- 



Cl. From (10.2) we find immediately 

S v (k) 



X 



f] 2 k 2 



(10.3) 



Thus the long range, Stokesian hydrodynamic propagation converts short range fluctua- 
tions in "P th into long range fluctuations in the fluid velocity. As mentioned in § 8J3 , the 
resulting divergence could lead to erratic coarsening rates and/or problems with finite 
size effects, throughout the viscous regime. This appears not to have been noticed by 
previous authors. 

There is a somewhat related anomaly that arises in colloidal suspensions under gravity, 
although in that case the short range fluctuations are in the density, which is effectively 
a random body force, rather than in a random stress: Segre, Herbolzheimer & Chaikin 



(1997|) 



10.2. Length scales from the velocity field 

The velocity structure factor, S v (k), can be used to calculate a velocity length scale, 
L V (T) analagous to L(T) (compare (4.1)): 

JS v (k,T)dk 



L V {T) = 2tt 



JkS v (k,T)dk' 



(10.4) 



This length measure was found to be insensitive to coarse-graining in nearly all cases. 
Data collected for L V (T) from the 256 3 runs was converted to reduced physical units, 
using the values of T int already obtained from the L(T) fits and given in table |[ (Hence 
no further fitting was involved.) The resulting scaling plot is shown, alongside the l(t) 
data presented earlier, in figure |ll] (left). The results in the inertial regime show a strong 
convergence between l v = L v /Lo and I = L/Lq: the velocity length l v shows the same 
i 2 / 3 scaling as I, with a similar prefactor. This is not obvious a priori, since, as mentioned 
above, the shapes of S v (k) and S(k) are very different. 
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More surprisingly, we find that to a fairly good numerical approximation, l v (t) shows 
a t 2 l z growth throughout the crossover region, and that this even extends far into the 
viscous regime, within which l v exceeds the domain scale I by a significant factor. How- 
ever, as the viscosity is increased to access the bottom left corner of the plot, the data 
are increasingly affected by finite size effects, since L v then is comparable to the system 
size A. These are especially pronounced for the most viscous run (with L v almost con- 
stant during that run). Allowing for these effects, the data is consistent with l v ~ t 2 ^ 3 at 
all times; however we have no reason to expect this result in the viscous regime, where 
both the simple and the extended scaling analyses (§ ||, [si) p redict instead l v ~ / ~ t. 



Note, though, that the velocity length measure chosen ( 10.4 ) is sensitive to the low k 
divergence found above in S v (k), which would give a contribution of order A/ In A from 
the lower limits of integration. It is possible that for the parameters and system sizes 
used here, this size-limited contribution combines with those from higher wavevector to 
give an apparent 2/3 power in the viscous regime. 

Alternative length measures may be had by taking the ratio of two other successive 



moments of S v (k), to replace ( 10.4 ). Adding one extra power of k to the top and bottom 
integrand gives a length measure that lies between l v and I throughout the viscous regime, 
reducing the exponent discrepancy there, but without attaining the linear scaling of I 
itself. For more than two extra powers of k the scaling gets worse, not better, as the 
integral in the denominator becomes dominated by high k contributions. 

10.3. Average velocities 

The RMS fluid velocities (spatially averaged) were calculated for all the runs in table 
H and are plotted in reduced physical units in figure [ll] (right) , alongside the reduced 
interface velocity l(t) derived from the order parameter. 

In the most viscous run, Run028, the RMS fluid velocity is larger than the interface 
velocity. Both velocities are fluctuating quite far from the expected constant behaviour 
in the linear region, and the fluctuations are more or less in step. This may in part be a 



facet of the erratic, finite-size limited behaviour seen in the far viscous regime (§ 3.5). 

Otherwise we observe that the RMS fluid velocity matches the interface velocity in the 
viscous and early crossover region, but grows larger than it in the inertial region, by 
about 40% at the largest l,t. The two most inertial runs (Run019 and Run032) appear 
to have the RMS velocity scaling with a slightly different exponent than the interface 
velocity (approximately as i -1 / 4 rather than i" 1 / 3 ), though this may not be significant. 
Such a deviation is foreseen by neither the simple nor the extended scaling theory, both 
of which have velocities scaling as v ~ I at all times. It would imply a buildup of kinetic 
energy in the fluid beyond that predicted by either scaling analysis. The excess may be 
caused by our approaching the limits of numerical accuracy in resolving velocity gradients 



with a consequent breakdown in energy conservation (see § 11.5 , 11.8 below). A similar 
breakdown, caused instead by having too high a viscosity in lattice units (as indicated 
in § 0), may likewise contribute to the excess RMS velocity seen in the most viscous run. 

10.4. Velocity distributions 

The PDF of the velocity components in a homogeneous, isotropic turbulent fluid is known 
to be almost Gaussian, and uncorrelated over large distances in both space and time. 



Non-Gaussian behaviour is found only in velocity increments and derivatives, (e.g. Monin 



k Yaglom 1975) 



The spinodal system is crucially different: it has a structural length scale L(T), and 
correlations and inhomogeneity can be expected at this scale. Hence the velocity com- 
ponents themselves can show non-Gaussian PDFs. The departures are expected to show 
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Figure 12. Left: flatness of velocity components for runs in table ^; runs with Ln > 0.5 (crosses) 
and with Ln < 0.5 (squares). Single fluid turbulence (filled triangles) are shown for comparison. 
Right: pdfs of velocity components for Run032 at time step 12000, with Gaussian PDF (solid) 
and transverse velocity derivative (dashed) shown for comparison. 



up through the fourth moment, since there is no preferred direction in space that would 
allow one to create a scalar from a third moment of velocity. (The third moment was 
checked and found to be close to zero.) The fourth moment is characterised by the flat- 
ness, (vt) / (Vn) 2 ~ 3, with a a Cartesian index; this vanishes for a Gaussian distribution. 

Figure [lj (left) shows the flatness for runs in table |^ as a function of time; a distinction 
is drawn between viscous runs (those with Lq < 0.5) and crossover/inertial runs (with 
Lo > 0.5). It can be seen that the flatness is quite variable, but as a general trend it 
grows slightly with time through each run, and also grows with decreasing Lq (increasing 
inertia). The velocity pdfs show correspondingly wider tails and narrower peaks than 
for a Gaussian; an example for Run032 is shown in figure |l2| (right). The shape is close 
to that found in the transverse velocity derivatives in the same system (shown dashed 
for comparison). These non-Gaussian effects are much more pronounced in the incrtial 
than in the viscous regime. 



11. Velocity derivatives 

We now turn to the analysis of spatial derivatives of the velocity field. Velocity deriva- 
tives come in two types, longitudinal, e.g., dv. x /dx, and transverse, e.g., dv x /dy or dv x /dz. 
Representative velocity derivatives (in fact, dv x /dx, dv v /dy, and dv y /dx) were calcu- 
lated. The differentiation was done by calculating ik x v y (\i) etc., and Fourier transforming 
back to real space. The derivatives are unambigous so long as the gradients of velocity 
are small on the scale of the lattice spacing, but otherwise the resulting gradient is not 



the same as taking a lattice derivative (which we define via (6.19) with v replacing 4>). In 



practice we have found signigicant (~ 40%) discrepancies the two methods; these were 



investigated in the context of energy conservation and are discussed further in § 11.5 For 
many purposes there is no particular reason to prefer the lattice to the Fourier definition 
and we retain the latter for simplicity. 

11.1. Skewness of velocity distribution 

In fully-developed turbulence the skewness of the longitudinal velocity derivatives ap- 
proaches —0.5, (e.g., Monin & Yaglom 1975). (The skewness of a variable y is (y 3 ) / (y 2 ) 3 ^ 2 .) 



The transverse derivatives have zero skewness, by symmetry, but positive flatness. 




Figure 13. Left: Skewness of the longitudinal velocity derivatives for Run031 (solid), Run015 
(dashed), Run018 (dot dashed) and single fluid decaying turbulence (long dashed). Right: Skew- 
ness of a longitudinal velocity derivative for Run031 (solid), a 96 3 run with the same parameters 
in which the interface was removed at time step 5000 (filled circles), and single fluid turbulence 
(dashed). The time scale for the turbulence data has been multiplied by 5 (left) or 10 (right) to 
facilitate comparison with the LB data. (The time dependence is not of relevance here once the 
initial stages are passed in both simulation methods.) 



Figure [13| (left) shows the skewness of the longitudinal velocity derivatives against time 
for the three most inertial runs in table |3[ Also shown for comparison is the skewness from 
the freely decaying turbulence simulation. In the two most inertial LB simulation runs, 
the skewness of the longitudinal velocity derivative reaches around —0.35. A plausible 
interpretation of this result is that patches of turbulence arise, but do not fill the whole 
system; if the patches have skewness —0.5, the overall value is less. From visualizations 
(see figure |l^) we know that the interface remains smooth so that any turbulent regions 
should be in the middle of the fluid domains. 

To test this idea further, a 96 3 run with the same parameters was done, and at timestep 
5000 (when the domain size was about 21, after residual diffusion had decayed to an 
acceptable level) the interface was suddenly removed by setting the order parameter to 
1 throughout the system. This converted it into a single fluid with the same velocity 
field, which was then allowed to evolve. The longitudinal velocity derivative skewness for 
this run is shown in figure [l^ (right), with Run031 and the single fluid turbulence data 
repeated for comparison. Once the interface was removed, the skewness quickly jumped 
to around —0.5 from —0.35, providing strong support for the 'turbulence in patches' 
hypothesis; it appears that on removal of the interface, the turbulence rapidly infects the 
whole system. (Fairly soon after this, the system became numerically unstable causing the 
skewness to rise rapidly back toward zero as spuriously large velocities were generated; 



this is visible in figure 13, but was not investigated further.) 



11.2. Reynolds numbers 
In this work we choose to regard all Reynolds numbers as estimates of the RMS ratio, 



i?2, of the nonlinear term to the viscous term in the NSE (3.2), as defined in (5.2). 

In practice, Reynolds numbers are usually constructed (via Re — TqIv/Lq) from a 
characteristic length (£) and a characteristic velocity (y). In a homogeneously turbulent 
single fluid, the usual choice of velocity scale is wrms- Also, the length measure i must 
itself be constructed from the velocity field; various choices then arise. One is the integral 
scale, defined (to within prefactor conventions) as L int = (3ir 2 p/2E) J kS v (k)dk where 



E = 2irp J k 2 S v (k)dk is the total kinetic energy. (This a close relative of our L v , (10.4.)) 
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Figure 14. Left: For runs in table [2| Re^, (bold line), and Re^ (diamonds) and for Runs in 
table §, Re A (dots). Right: Ratios Ri and R 2 for runs with L = 36, 2.9, 0.59, 0.15, 0.054, 0.024, 
0.01, 0.01 (different parameters), 0.0016, 0.00095, 0.00039, 0.0003. System sizes are 96 3 (open 
symbols) and 128 3 (filled symbols). 



A second is the Taylor microscale, A, denned in Eq. (5.7). A third length is the Kol- 
mogorov microscale, A^, ( |5.§| ), but this is not normally used to form a Reynolds number. 
In practice, most isotropic homogeneous turbulence simulation studies use as Reynolds 
number Re\ — pvbms^/v- This should give a reasonable estimate of the ratio i?2, be- 
cause as noted in § ^, A is the length scale associated with the V operator in Vv. But 
note that, according to the extended scaling analysis (table Q), it is still not an accurate 
estimate asymptotically at late times: there one has the prediction that Re\ ~ Aw ~ T 1 / 6 
while R 2 ~T°. 

In the spinodal system an obvious alternative to Re\ arises from choosing L(T), the 
domain size, as the length scale; either urms or the interface velocity L can then be used 
for the velocity scale. However, the resultings Reynolds numbers, Re^ = pifRMs/ 7 ?! an d 



i?e</> = pLL/i] (see (5.4)) are not directly comparable with Re\, which is always much 
smaller. 

All three quantities are shown as functions of the reduced time t in figure [l4| (left). 
There is little difference between Re<j> and RtL\ both show reasonable scaling, with the 



deviations in the most extreme runs stemming from those discussed already in § |10.3 



The large t asymptote for Re^ in the inertial region is approximately Re^ ~ i 1 / 3 as 
predicted by both the simple and the extended scaling analysis (but questioned by prant 



fc Elder 1999| ). Our simulations extend from 0.1 < Re^ < 350, and the crossover region 
occupies the range 1 < Re^ < 100. Thus in terms of Re^ (rather than t) the crossover is 
not, after all, quite so broad. 

In contrast to Re^ and Re^, the data for Re\ does not show good scaling behaviour. 
The overall trend is of linear scaling in the viscous region and slower growth at around 
i 1 / 6 in the inertial region; both are broadly consistent with the extended scaling theory 
(see table Q). However, the individual runs do not line up onto a single curve. This non- 
scaling behaviour of Re\ can be traced to that of A itself, which is examined in more 
detail in § |Tf , 

11.3. Ratios of terms in the NSE 
Shown in figure [l4| (right) are the actual RMS ratios Ri and R2 defined in (5.1) and 



( p.2\) respectively. (Recall that the latter is the ratio of nonlinear to viscous terms in the 
NSE, which is what we believe a Reynolds number should estimate.) In order to form 
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these ratios, one must first evaluate (RMS values of) the three relevant NSE terms, ?yV 2 v, 
pv and p(v.V)v, which are vector quantities at each point in space. An RMS value is 
calculated as, |/|rms = ({fx) + (/«) + (/z)) 1/2 ; where f x , f y , and f z are the Cartesian 
components of the vector f = t)V v (for example), and the average is taken over the 
whole system. The first and second order spatial derivates of the velocity were found 
from fc-space data (as in § [TTl] ) on 96 3 and 128 2 runs; a further complication is that to 
evaluate v, velocity data from consecutive time steps is required. Due to data storage 
constraints, this was only collected for a set of 96 3 runs. This is large enough to obtain Rx 
values for a domain size just larger than L m i n , but not to determine its time dependence 
accurately within any particular run; checks were made for consistency by comparing 
other quantities with 128 3 data. 

The values of R\ and R2 shown in figure |IJ (right) lie between Re^ and Re\ but are 
closer to the latter. The two ratios remain the same order of magnitude as each other 
throughout, but vary by three orders of magnitude, from R2 ~ 10~ 2 at the viscous end 
(indicating the viscous term is dominant by two orders of magnitude) to R2 ~ 10 1 at the 
inertial end (indicating that the inertial terms are dominant by one order of magnitude) . 
Though we cannot calculate it directly, R2 is presumably somewhat higher (about 20) 
towards the end of our most inertial 256 3 run. This confirms the claim made in § ^| 
that our simulations have reached a regime where the inertial terms are dominant in the 
dynamics. 

Looking more closely at the behaviour of R\ relative to R2, there is a significant 
difference in the crossover region, by around a factor of two (i?i > ffe)- Then, in the 
inertial regime, R\ becomes less than R2 by about 50% and appears to be heading for 
a lower growth rate. This deviation suggests that the asymptotic behaviour of these two 
ratios is perhaps going to be different. That would be consistent with the extended scaling 
predictions which are (table |l|) that R\ ~ t~ Y ^ while R2 — > constant; but if so neither 
curve is close to its final asymptote yet. 

That the LB simulations are still far from the final asymptotic behaviour in our most 
inertial runs, even though inertial terms are clearly dominant there, is not unreasonable. 
The largest ratio R2 — 20 actually achieved in our simulations is, by tur bulence standar ds, 
still a fairly low Reynolds number. If the extended scaling theory of Kendon (2000| ) is 



correct, the final asymptotes for quantities relating to fluid velocity cannot be attained 
until an appreciable 'inertial range' has developed between the interfacial driving scale 
L(T), and the smaller scales (A, Xd) where energy dissipation is taking place. However, 
to observe scaling of l{f) for the structural (as opposed to velocity) data, it may not 
be necessary that the inertial cascade is fully established; one might only require that a 
reasonable degree of decoupling between interfacial motion and viscous dissipation has 
taken place. Our results for 1(f) suggest that this has already happened by the time 
R 2 = 10. 

11.4. Structure factors of the NSE terms 

Further information on the behaviour of the NSE terms can be obtained by calculating 
the structure factor for each term. (These are (|f(k)| 2 ) where f = ?yV 2 v,v.Vv or pir.) 
Results for one viscous and one inertial run are shown in figure |l5[ Looking first at 
the viscous run, the structure factor of the viscous term takes the form of a broad 
peak stretching from a small bump at wavelengths around 12 lattice units (which is 
around half the domain size, L(T) = 25), down to the interface width, £ ~ 3. Thus the 
dissipation is taking place over the smaller length scales in the system. (The small bump 
is a manifestation of the domain size in the dynamics; it is present at around L(T)/2 
throughout the run.) 
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Figure 15. Structure factor of the viscous (squares) and nonlinear (triangles) terms in the NSE. 
Left: Run027 is in the viscous regime with Ln = 36. Right: Run031 is in the inertial regime with 
Ln = 0.0003. Both are shown on a log-log plot for the timestep at which L ~ 25 lattice units. 
L — 25 is marked on the x-axis along with the interface width, £ ~ 3. (The pdv/dT data is 
from the corresponding 96 s runs.) 



In the inertial regime, the viscous term is, as expected, smaller than the inertial terms. 
The viscous term is similar in shape to that in the viscous regime, with the addition 
of two large peaks at high k that presumably arise from the presence of the interface 



(compare § |8.4| ) and/or lattice effects. This suggests perhaps that the largest dissipative 
forces, and therefore most dissipation, are happening close to the interface (or at least 
on that length scale). The acceleration term has a broad peak in the structure factor 
at wavelengths around L = 25, and tails off quite sharply below 10 lattice units. The 
nonlinear term has a structure factor with a broad peak centred around 15 lattice units, 
intermediate between the length scale of the interface, L(T) = 25, and the dissipation 
length scales. 

The overall picture, though not quantitative (in view of the numerous sources of uncer- 
tainty, especially in the velocity derivatives) is qualitatively consistent with the extended 
scaling picture, in which the nonlinear term transports energy from large to small scales, 
where it is then dissipated. But in these runs there is still considerable overlap between 
the length scales for each term, as expected at relatively low Reynolds number. 

11.5. Energy conservation, dissipation rate 

The dissipation rate is crucial to the energy balance in the simulation. The LB algorithm 
is isothermal, and therefore does not strictly conserve energy, in that the energy dissipated 
as heat by the viscous stresses is taken out of the system locally rather than having to 
diffuse thermally to the boundaries. Nonetheless, the simulation should faithfully recreate 
the energetics of the isothermal Navier-Stokes equation, driven by interfacial motion (as 



in (3.2)), which makes precisely the same assumption. In the inertial regime, the inter- 



facial energy should first be transformed into kinetic energy of the fluid; the observation 



of I ~ t 2 / 3 (see § EL3) shows that this is happening at roughly the expected rate. To 
complete the energy balance, energy must also be removed from the simulation at the 
correct rate through viscous dissipation. This requires accurate modelling of the smaller 
length scales (A — Xd) at which such dissipation is taking place. 

The dissipation rate can be calculated directly from the velocity data as 

e = )j((V a t) ? )(V a ^)) = r}f k 2 S v (k)d 3 k (11.1) 
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Figure 16. Energy conservation in viscous region. Sin (circles), e (squares), for Run008 (filled), 
Run014 (open), and Run027 (solid line). In these runs v is negligible on the scale of the graph. 



where the angle brackets are spatial averages and both the incompressibility (V a v a = 
0) and homogeneity ({v a vp) — (l/'3)S a pS v (k)) of the binary fluid system have been 
exploited. The second expression is tantamount to using the Fourier definition of spatial 
derivatives, and as mentioned previously, this need not coincide with lattice derivatives 
when gradients are not small. Therefore the dissipation rates were calculated first using 
the Fourier space expression in ( |ll.l[ ) and then, for some but not all runs, using the real 
space expression with the lattice gradient operator as defined in Eq. ( 6.19 ). 

The two estimates (Fourier and lattice) of the dissipation rate were compared with 
the rate of decrease of the sum of the interfacial energy and the fluid kinetic energy 
in the system: these should balance, once the interfaces are well formed and residual 
diffusion has become negligible. (Note that diffusion introduces its own contribution to 
the dissipation, not accounted for in the NSE, which may dominate early in each run.) 
We found that the lattice derivatives gave better agreement than Fourier ones in this 
comparison. However, even with these, the results were never more than satisfactory; 
best in the crossover region (within 25% of the expected value) but giving dissipation 
rates well above those required by energy conservation in the viscous regime, somewhat 
below in the inertial. At the viscous extreme the discrepancies were a factor three or 
more. Even in the crossover region, there was a tendency for the computed dissipation 
rate to drift below that required for energy conservation. Examples are shown in figure [l6[ 
Use of Fourier derivatives to estimate the dissipation rate instead gave values that were 
consistently too high, through all regimes, by a factor 2 or so for runs in the crossover 
and inertial regimes and much more than this in the viscous regime. 

A check was made to ensure that these discrepancies in the dissipation rate did not arise 
from compressibility effects, by computing the full expression for e, without assuming 
fluid incompressibility as was done in (11.1). (See Landau fc Lifshitz (1959j ) for the 
relevant expression.) This check was done with lattice derivatives; the results did not 
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Figure 17. Left: Dissipation rate at time step when domain size is L = 30 (in simulation units) 
for Runs in table || (circles, solid line) and Runs in table (plus, dotted). Right: Dissipation 
scale, Ad, for runs in table H (plus, dashed), and runs in table tfl (stars, solid); Taylor microscale, 
A for same runs (diamond, dashed; circles, solid respectively). All in lattice units with domain 
size L = 30 lattice units. 



differ appreciably from the lattice version of Eq. (11.1), so that compressibility is not 
responsible here. 

The fact that the two methods of finding derivatives differ, shows that dissipation is 
primarily taking place at short scales (of order a few lattice spacings). The results do 
not imply that LB is failing to conserve energy, but do show that the viscous dissipation 
actually generated by the algorithm is not accurately estimated from the Fourier deriva- 
tives; the lattice difference estimates for e are better, but still not accurate. In principle 
even these need not give the true dissipation since the LB algorithm actually dissipates 
energy by relaxing the velocity distribution functions ( |6.5| ) , and not by calculating lattice 
gradients. Note, in any case, that for qualitiative comparisons between runs (requiring 
log-log plots covering many decades) the factor two difference between Fourier and lattice 
derivatives is barely detectable, and we ignore it, for those purposes, below. 

For the dissipation rate itself, simple scaling theory predicts that e will always scale 
as t~ 2 . Kendon (200C| ) predicts the scaling to be e ~ t~ 2 in the viscous region, but 
e ~ £-5/3 j n tfic incrtial region (table [I]). A scaling plot of the (Fourier) dissipation rate 
(in reduced physical units, as usual) is shown in figure [ij] (left). To the accuracy we are 
working, these data barely discriminate betw een the two sca ling predictions though the 
upward curvature may slightly favour that of Kendon (2000| ). There is some sign in the 



viscous regime that the rate increases too fast; this might be connected with the velocity 



anomalies discused previously (§ 10.1) 



11.6. Taylor and Kolmogorov microscales 



The Taylor and Kolmogorov scales were defined in (5/7) and (5.8) respectively. Figure 
|l7| (right) shows A, Xd in simulation units for each run in table |2| and table [| Values are 
evaluated mid-run at time T such that the domain size L(T) = 30, but plotted against 
the reduced physical time. The curves thus have the same shape as a scaling plot of 
\/L(T) and Xd/L(T) (since L(T) is here fixed in simulation units) but also allows the 
A, Ad values to be compared with the lattice spacing (unity) and the linear system size 
(128 or 256). Given the expected effects of coarse-graining on the 256 3 data at short 
length scales, the agreement between the two lattice sizes is satisfactory but only the 
128 3 data is discussed in what follows. 
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As can be seen in figure [l?] (right), X/L(T) does not vary by much between runs, al- 
though it is roughly a factor of two larger in the crossover region than in the viscous and 
inertial extremes. The extended scaling analysis (table |l|) predicts that X/L(T) ~ i -1 / 6 
in the inertial region, and the data is broadly consistent with this; an asymptotically con- 
stant ratio, predicted by simple scaling, is marginally less plausible though not ruled out. 
Likewise, X c i/L(T) is not far from the prediction of the new scaling theory (i -1 / 4 ), and 
perhaps somewhat further from that of simple scaling (t -1 / 6 ). For the viscous regime, 
however, Xd/L(T) does not appear to scale as i -1 / 2 , nor does X/L(T) approach a con- 
stant, as predicted by both scaling theories. Once again the low k velocity anomaly (§ 



10.1 ) in the viscous regime c ould be to blame for this, but so could inaccuracies in the 
estimated dissipation rate (§ 11.5| ). 



11.7. Resolution of the energy cascade 

The Kolmogorov microscale has a further significance for any attempt at numerical sim- 
ulation in inertial fluids. In a fully turbulent fluid, Xd is expected to be smaller than A, 
and to mark the smallest length scale relevant for dissipation. This should, if the results 
are to be relied upon, be 'resolved', that is, Xd should lie above the discretization scale 
set by the lattice. There is some debate over exactly how small Xd can be in relation 
to the lattice spacing, but a factor of 1.0 to 1.6 h as been put forward; see for example, 
Eswaran fc Pope (1988j ); Yeung fc Brasseur (1991 ). Our best estimate for Xd in the most 



inertial (256 3 ) run is around A^ = 1.8 lattice units; values are higher in all other runs. 
H owev er, in view of the quantitative uncertainty in our velocity derivative estimates (see 



§ 11.5 ), this estimate is good to only a factor of 2 or so. Therefore, the most inertial of 
our 256 3 runs is certainly close to the resolution limit, and we would not wish to proceed 
to higher l,t without a larger system size. Note that only in that run is Xd appreciably 
smaller than A as expected asymptotically in a turbulent flow; but the Taylor scale A is 
well-resolved, by the same criterion, in all runs. 

However, estimates of how well the cascade is resolved using these two lengths could be 
misleading, in that they assume local as well as global homogeneity. So if, for example, 
most of the dissipation were to occur in a thin layer around the interfaces in the system, 
the globally averaged Xd might suggest that the dissipation was fully resolved whereas 
in fact it was not. This would be consistent with our findings (§ 11. 5| ) that the actual 
dissipation rate in LB is poorly estimated from the Fourier velocity derivatives, and 
imperfectly even from lattice ones. 

On balance, we believe that the energy cascade is adequately resolved in most of 
our LB simulations, at least for the purposes of getting the correct evolution of the 



structural length scale 1(f). (As noted in § 11.3, a fully resolved cascade might not be 



needed for this, so long as there is sufficient decoupling of interfacial motion and viscous 
dissipation.) However, the most inertial runs may well be marginal in terms of resolution; 
in common with several other aspects of the simulations (such as the residial diffusion, 
and anisotropy), these runs are at the limit of what is possible using this simulation 
method with current computational resources. 



11.8. Apparent scaling violations 

In our various analyses of the velocity derivatives, we focused mainly on scaling plots 
made by taking a representative 'mid-run' data point from each LB simulation run (e.g. 
figure [It], left). These plausibly connect to form a smooth curve. However, one test of 
scaling is whether the full datasets (rather than a representative point) from different runs 
appear to join up smoothly when plotted in reduced physical units; this was satisfactorily 
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t = 0" - T (nt ) / T t = (T-T lnt )/T 



Figure 18. Left: NSE terms for all runs scaled. Most of this data is for 96 3 systems, with 
128 3 data where available. Right: Dissipation scale, Ad, for runs in table ra (open circles), Taylor 
microscale, A for same runs (shorter solid lines) and runs in table |^ (filledcircles) . All in reduced 
units. 



the case for the structural data l(t) (figure ||) and reasonable also (though not by any 
means perfect) for the velocity data l v (t) (figure |TT|, left), Rcl and Re^, (figure [lj], left). 

However, the expected scaling within runs failed completely for Re\ (also in figure [l4], 
left). This quantity differs from the others just mentioned in that its definition involves A, 
which in turn involves the dissipation rate e (see (5.7)) and therefore depends on velocity 
derivatives. More generally, we find similar 'scaling violations' in all the quantities we 
have looked at that involve calculating velocity derivatives. Figure fl8| shows this for (left) 
the various RMS terms in the NSE and (right) the Taylor and Kolmogorov scales. The 
dissipation rate (not shown) shows similar features. 

Apart from the failure of the runs to join up, the data for the NSE terms is in reasonable 
accord with the scaling predictions (table |l|) for the viscous term of t~ 2 in the linear region 
and <~ 5 / 3 or i~ 7 / 6 in the inertial region, and for the inertial terms of i~ 4 / 3 or £~ 7 / 6 in 
the inertial region. Even the scaling of t^ 1 for the nonlinear term in the viscous region 
is as predicted by the new scaling theory (Kendon 2000). However, within each run the 
quantities are falling more slowly than all these predictions. Similarly A and Xd do not 
scale within single runs in the same way as they scale between runs, except perhaps in 
the middle of the crossover region. (The same is true of e, the dissipation rate, from 
which Xd is derived.) 

At present, we have no simple explanation of these apparent scaling violations. They 
could perhaps be a sign of subtle nonuniversalities of the type suggested by [Jury et ~al 



(19996| ), though it would certainly be premature to conclude this without similar (though 
perhaps less severe) effects being detected in l(t). Another possibility, noted previously (§ 



11.5, 11.7) is that the we are not able, using LB with these system sizes, to fully resolve 



the velocity gradients that arise. 

On the other hand, our velocity derivatives do not appear to be completely lattice 
controlled, since for example, in (11.1) one could then set Vv ~ v. This would give 
£ *~ L , giving in the inertial regime X d ~ T 1 / 6 within each run (which is about right) 
but in the viscous regime A^ ~ T°, and also A ~ T° throughout both regimes, neither 
of which is plausible for our data (see figure [l8], right). We have made various further 
attempts to explain the data by assuming 'composite' scalings, based for example on the 
idea that dissipation occurs mainly near the interface where velocity gradients might be 
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anomalously large (Vv ~ i/OR ^ n that case, one would have e ~ (£/L)i 2 /£ 2 , where 
the factor £/L is the volume fraction of interface. Although none of these attempts was 
successful, an explanation along these lines is not completely ruled out. 

One further complicating factor is that, within LB, the fluid viscosity is frequency 
dependent: the Navier Stokes limit requires that velocities are not only small (in lattice 
units) but also evolving sufficiently slowly. For very low viscosities, t% approaches 1/2 
which means that the distribution function fa is over-relaxed (nearly reversing sign at 
each time step), see Equation 6.6. (In effect, the 'natural' viscosity scale in LB is of order 
unity and the over-relaxation is used to create an unnaturally small value.) This works 
successfully if accelerations are small but could result in a higher dissipation rate than 
would be the case for a purely newtonian fluid, especially in the most inertial runs. 



12. Conclusions 

Our LB results for symmetric binary fluids, undergoing spinodal decomposition after 
a deep quench, appear to confirm the dynamical scaling hypothesis, which requires a 
universal time dependence of the structural length scale L when expressed in reduced 
physical units as a function of time, l(t) (figure |9|). We found no signs of nonuniversality 



here, although some weak breakdown of it, as suggested by Jury et al. (1999&), cannot 
be entirely ruled out. By exploitation of the expected scaling, and careful parameter 
steering and validation tests to eliminate residual diffusion and other unwanted effects, 
we achieved an l(t) curve spanning seven decades of reduced time t and five of reduced 
length I. The l(t) curve asymptotes to I = bit at t <C t* (the viscous regime), with 
bi ~ 0.072, a nd I = 6 2 /3^ 2/3 at t > t* (the inertial regime), with b 2 / 3 ^ 1.0. On 



our reading of Gucnoun et al. (1987 ), b\ is within 10% of the most careful expe rimental 
measurements (albeit not for a deep quench) although others ( Laradji et al. 1996] ) extract 



an estimate about twice as large from the same experimental results. 

The crossover time, t* in reduced physical units, as defined by the interception of the 
viscous and inertial asymptotes, is surprisingly large (t* ~ 10 4 ). So is the width of the 
crossover region (four decades). This may explain why the t 2 ^ region has never yet been 
confirmed in laboratory experiments. The default assumption that t* is 'of order unity' 
would mean that for fluid pairs typically studied (short chain alcohols and water for 
example), the inertial regime could be accessed at length scales of a few microns, easily 
studied by light scattering. But in fact, to exit the crossover region (at t* ~ 10 6 ) one 
needs I = L/Lq > 10 4 which for typical fluids gives L of order centimetres (and larger still 
for the near-critical quenches often used, see Gucnoun et al. (1987[ )). This requires direct 



visualisation experiments, not light scattering, and more importantly requires rigorous 
exclusion of thermal convection and gravitational effects. The latter is likely to be possible 



only with extremely careful density matching, or in microgravity, Cambon et al. (1997). 

Our findings are less extreme when translated into Reynolds numbers. The conventional 
definition of the Reynolds numbers for the spinodal problem is Re^ — 11; the crossover 
regime then spans 1 < Re^ < 100 and the Reynolds number range we actually achieved 
was 0.1 < Re^ < 350. (Note that we see no sign of a saturating Re^ as predicted by 



Grant & Elder (1999), but cannot rule this out, at some value well above 350. ) 



f We note that in a very recent analysis, Solis & de la Cruz (2000) developed an heuristic 
alternative to that of Kendon (200Cp in which dissipation takes place in an (asymptotically) thin 
layer around the interlace rather than in the bulk. They argue that coarsening is limited by the 
damping rate of capillary waves of wavelength L, and that this damping rate is unaffected by 
nonlinearity; this gives I ~ t 4 ' 7 . We find no evidence for this in our data. Using the oscillation 
frequency instead of the linear damping rate recovers I ~ t 2//3 . 
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These figures are reduced further if the actual ratio R2 of nonlinear to viscous terms 
is used as the Reynolds number; our results span 0.01 < R2 < 20. At the upper end of 
this range, we have described in detail the first unambiguous simulations of a regime in 
which the inertial terms in the Navier-Stokes equation are actually large compared to the 
viscous ones; here we observed I ~ i 2 / 3 . However these simulations may still be far from 
the asymptotic regime addressed by the scaling analysis of Kendon (2000] ) , in which ulti- 
mately one expects <C A <C L for the Kolmogorov, Taylor and structural length scales. 
But in fact, although A^ is only slightly smaller than A even for our most inertial run, our 
data for the inertial regime is already somewhat more consistent with Kcndon's a nalysis 
than with the simpler (single length scale) scaling predictions of Furukawa (1985 ). 

In the viscous regime, we recovered the expected (I ~ t) scaling though this is somewhat 
erratic and accompanied by an unexpected behaviour of the velocity correlation length 
l v = L v / Lq (figure |]J left) and possibly also the dissipation rate (figure 17, left). The first 
of these, at least, may be related to the presence of anomalous long-range correlations 
of the velocity over length scales much larger than the domain size L(T). These were 
detected numerically in our most viscous runs, and can be explained simply by assuming 
that the interfaces contribute a random stress with local correlations on scale L (see § 



10.1). It is quite possible that such effects have arisen undetected in previous simulations 



of the viscous regime; in principle they could lead to finite size problems arising long 
before the domain size L approaches the size A of the simulation cell. Fortunately, the 
same does not appear to happen in the inertial regime. 

As well as through their differing domain growth laws (I ~ t or I ~ i 2 / 3 ) the viscous 
and inertial regimes can be distinguished by rather subtle changes in interface geometry 
(evident through the structure factor and by visualisation, with more thin necks in the 
inertial case) and by much larger changes in the velocity statistics. For example, the 
velocity PDF is significantly flatter for inertial than for viscous runs and also there is 
significant skewness for longitudinal velocity derivatives in the inertial regime. In that 
regime, we achieve results that show some of the characteristics of a turbulent fluid, 
although non-Gaussian features are seen directly in the velocity distribution as well in 
that of derivatives. This reflects the presence of microstructure (interfaces), as well as 
some turbulence, in the fluid mixture. 

The scaling of quantities involving velocity derivatives was found to give reasonably 
continuous curves when mid-run values were plotted in reduced physical units (e.g. figure 
|l7| , left) but showed apparently systematic violations of the scaling behaviour which that 
would imply, when analysed in detail within each run. This could indicate some physical 
nonuniversality, or some shortcoming of the LB code; but equally it can be attributed to 
uncertainties over how to accurately estimate derivatives that are not small on a lattice 
scale. The LB code is implicitly aware of differences in the local velocity distributions 
but does not itself calculate velocity derivatives. Therefore, when these derivatives are 
not small, there can be discrepancies between their 'effective' values (as defined, e.g., 
through the actual dissipation rate generated by the algorithm) and any estimate based 
on Fourier, or lattice difference, data. 

The fact that such ambiguities arise at all is an indication that our LB simulations are 
close to the limits of accuracy acceptable for the fluid motion, since ideally all gradients 
encountered in LB should be small, not large, on the lattice scale. It is interesting that 
they arise even when the Taylor and Kolmogorov microscales, A and Ad, are both signifi- 
cantly larger the lattice spacing; this perhaps suggests heterogeneity of the velocity field 
that is stronger than in single-fluid homogeneous turbulence. It may mean that, despite 
our best endeavours, we have not adequately resolved the fluid motion at short length 
scales. Fortunately, according to the analysis of Kendon (2000| ), the main requirement 
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for observation of I ~ t 2 ^ 3 in the inertial regime is that the interface motion is decou- 
pled (by fluid convection) from the dissipation scale, and not that the latter is modelled 
with complete accuracy. The reason is that, once the interface is decoupled, its motion 
is controlled by how fast it feeds kinetic energy into the fluid at large scales; the pre- 
cise details of the dissipation further down the cascade has no further influence on the 
interfacial dynamics. If so, statistics based on the domain growth law l(t) may be much 
more robust than those for the fluid velocity and (especially) velocity derivatives, in the 
inertial regime. 

Very similar accuracy limits were also reached in the thermodynamic sector, where we 
required narrow interfaces, with large composition gradients, to satisfy the conflicting 
requirements of rapid interfacial equilibration and low residual diffusion. (This gave in- 
terfacial tensions f 0-f 5% different from nominal, with measurable lattice anisotropy.) We 
have taken unusual care to analyse (and maintain under reasonable control) the various 
sources of error. The reader has, we hope, enough information to decide for herself how 
much confidence to place in our various results. 

It is interesting to ask how simulations might be taken further into the inertial regime. 



If the analysis of Kendon (200C ) is any guide, the ultimate asymptotic state entails, 



computationally, resolution of the following scale hierarchy: 

a<£<A d <A<£<A (12. 1) 

with a the lattice spacing, or, more generally, a -3 the density of degrees of freedom in 
the simulation. A factor 3 between each of these length scales is roughly what we achieve 
at A — 256. But a worthwhile improvement to give, say, a factor 10 between each would 
require A = 10 5 . For a three dimensional system, this lies beyond any foreseeable inno- 
vation in computational hardware if any current methodology is used. Possibly the way 
forward would be to couple a coarse-grained algorithm (such as large eddy simulation) 
for the turbulent fluid to a moving interface; we leave this to others to explore. 

We would like to thank Alan Bray, Anatoly Malevanets, Alexander Wagner, Patrick 
Warren, Julia Yeomans and Alistair Young for helpful discussions. Work funded in part 
under EPSRC GR/M56234. 



Appendix A. Analysis of hydrodynamic modes 

In this appendix we analyse the hydrodynamic modes for the LB hydrodynamic equa- 
tions (6.21-6.22). We take as our reference state one in which the fluid is at rest with 
uniform concentration. In this case, to linear order in the perturbation of the hydrody- 
namic fields the Navier-Stokes equation is satisfied, and only the spurious terms in the 
balance equation for the order parameter remain. Specifically, the equations we are going 
to analyse are 



V.(pv) = 



v = — V/P 
P 

V.(0v) - (r 2 - ~) 



th 



iA7 2 v + CV(V.v) 



MVV-V.f-V.?" 1 ) 
P 



(Al) 



where v = r\j p and £ = (£ + 2rj/3)/p are the kinematic shear and bulk viscosities, 
respectively. Using eq.(2.9), the last term in the third equation can be rewritten as 



44 V. M. Kendon, M. E. Cates, J-C. Desplat, I. Pagonabarraga and P. Bladon 



V .V th = c 2 s V p + <j)V p (A 2) 

where c s is the speed of sound and p the chemical potential. This expression shows 
that there is a spurious coupling between the order parameter and gradients of the 
density field. Such a coupling can only be relevant when the binary mixture exhibits 
compressibility. In terms of the deviations of the hydrodynamic fields with respect to 



their equilibrium values, po = 1, </>q, v = 0, eqs.( A 1) become 
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where we have defined u>2 = T2 — 5- In Fourier space, this set of equations can be expressed 
in matrix form, setting X = (Sp, k.v, 6(f), <5v — kk.<5v/fc 2 ) as 
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where a = A + 3B(f>Q. The hydrodynamic mode eigenvalues are correspondingly 
Ai = -vk 2 
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The first mode corresponds to the propagation of shear waves, the subsequent two 
modes are related to the propagation and damping of compressible momentum waves, 
and the final mode is related to the diffusion of the order parameter. One can see that the 
presence of the spurious term in the convection-diffusion equation does not qualitatively 
modify the propagation of the order parameter: it remains diffusive. However, its eigen- 
value depends on the compressibility of the the fluid, and for a compressible fluid there is 
then a reduction of diffusivity coming from the coupling to sound waves. Such a coupling 
also modifies the propagation and damping of longitudinal waves in the medium, so that 
these become a function of the order parameter. 

To whatever extent the compressibility of the fluid can be neglected, the additional 
coupling between the order parameter and the density becomes irrelevant, and the ex- 
pected diffusion coefficient for the convection-diffusion equation of the order parameter 
is recovered. 
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Figure 19. Scaling plot in reduced variables (L/Lo, T /To) Bold lines (left to right) are LB 256 3 
data from table (t op to bott om). Also shown are results fr om other published work: squares 
Appcrt et al. (199E), triangles Laradji et al. (199£| ), circles Bastea & Lebowitz (1997), Inset: 
DPD data of Jury et al. (19996) (solid lines) with one LB data set (Lo = 0.15, pluses) repeated 
for comparison. 



Appendix B. Comparisons with other work 

In order to compare our results with others' published work, a similar scaling procedure 
must be applied to place their data onto the universal scaling plot in figure ^|. This 
is only possible for work reported in sufficient detail to enable values for Lq and To 
to be calculated. Recent three-dimensional work where comparison is possible includes 
that of [Bastea fc Lebowitz (1997 ), Laradji et al. (1996), |Appert et al. ( 1995| ) , and lury 
et al. (19996). These four sets of results are shown along with the LB data in figure [19]? 
Simulations of three-dimensional spinodal decomposition with hydrodynamics for which 



quantitative comparisons were not possible include, Koga & Kawasaki (1991), Puri & 



Diinweg (1992), and Alexander et al. (1993), all of whom claimed to have simulated the 
linear regime; Shinozaki fc Oono (1991[ ), [Ma et al. (1992| ), |Lookman et al. (1996| ). The 
last two claimed to have simulated the inertial regime but offered no evidence beyond 
their fitted exponent values for definitely having inertial rather than diffusive effects. 

The four studies for which detailed numerical comparisons are possible will now be 
considered in turn. 
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B.l. Bastea & Lebowitz (1997) 



Bastea & Lebowitz (TsW) carried out a three-dimensional simulation containing about 
1.4 x 10 6 particles whose motion was described mesoscopically by Boltzmann-Vlasov 
equations. They combined direct simulation Monte Carlo methods for the short range 
interaction with particle- in-cell methods for long range interactions. The fluid system 
is relatively low-density; they describe it as a gas-gas phase sep aration. The necessary 
fitting and scaling for the results reported by Bastea fc Lebowit^ was done by Jury et al 



19996|) . 



On the universal scaling plo t, figure [19|, the supposedly viscous-regime portion of the 
data from Bastea <fc Lebowitz is shown as circles. It lies in the correct (lower left) region 
of the graph, but well to the left of any of the LB results. The scaled value of the fit 
parameter, b\ = v/ (Lq/Tq), from the refitting by Jury et al. is b± = 0.3. This compares 



with the considerably smaller values from the LB results in the linear region. The dy- 
namical scaling hypothesis requires however that b\ should be universal to all systems in 
the viscous regime. The high reported b\ leads us to suspect residual diffusion. Bastea 
fc Lebowitz do not repor t th e diffusion rate in their system in a form that can be used 
to apply the analysis of § 8.3, but test LB runs have been done with high diffusion rates 



that produce LB data sets very similar to those of Bastea & Lebowitz. Two such test 
runs, Run024 and RunOlO, are shown in figure [20| 

On the left, Run024 is compared with Run028, one of the runs in table ||. Run024 
has the same parameters as Run028 except for the mobility, which is eight times higher. 
Run024 fits a free exponent of a = 0.85, which is similar to Run022 (also in the linear 
region). However, the lack of an initial flat diffusive region (also absent in Bastea & 



Lebowitz (1997)) strongly suggests the diffusion is too strong for uncontaminated linear 



growth to be observed within the attainable system size. A linear fit to the upper part 
of the data produces a scaled value of b\ — 0.082, compared with 0.073 for Run028. 

On the right, RunOlO has a value of Lq = 381, larger than for any of the runs finally 
used by us, see table ||. This should put it even further into the linear region than the 
rest of the runs. However, a fit with a free exponent produces a = 0.7. A linear fit to 
the upper part of the data produces a value for the scaled fit parameter of b\ = 0.32, 
i.e. about the same as the data from Bastea & Lebowitz, which this curve resembles. It 
seems likely, therefore, that the data from Bastea & Lebowitz (1997) has strong residual 
diffusion, and the results they present are a mixture of linear and diffusive growth. 
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B.2. 



Laradji, Toxvaerd & Mouritsen (199t) 



Laradji et al. (1996) used a large-scale molecular dynamics simulation of a Lennard-Jones 
model with 343000 particles with a deep quench. The necessary fitting and scaling for 



the data of Laradji et al. was done by Jury et al. (19996). The results are shown on 
figure 19 as triangles lying, like those of Bastea & Lebowitz (1997), to the left of the 



LB data towards the lower left corner, again in the viscous regime according to the LB 



analysis. Laradji et al. claimed their results confirmed the linear scaling, but their value 
of the scaled fit parameter, b\ — 0.13, is again higher than ours. However, it is in well 



within the range of values spanned by the the over-diffusive LB runs shown in figure 20 
(&i = 0.082,0.32). The shape of their L(T) curve matches that of Run024 in figure |0 
(left), i.e. there is no initial diffusive plateau before domain growth begins. In LB data 
this was always a sign that the run had too high a mobility and significant contamination 
by residual diffusion. 



B.3. Appert, Olson, Rothman & Zaleski (199l\ ) 



Appert et al. (lyyb) used a three-dimensional lattice gas simulation to simulate spinodal 
decomposition. Their largest system size was 128 3 and they also used 64 3 to test for 
finite size effects; they rejected data with L > A/2, whereas we saw significant finite size 
effects by this stage (and applied the stricter condition, L < A/4). They claimed a fitted 
exponent of a ~ 2/3, which, if correct, would put their results in the inertial regime. 
Taking the relevant interval of their data, refitting it by our method (giving an exponent 
a = 0.62) and converting into reduced physical units gives the dataset (s quares) on figure 
fL9[ This shows that, if our LB data is correct, the data of Appert et al. is actually in the 
crossover region. It appears to be asymptoting onto the LB data from above, which again 
suggests significant residual diffusion, giving too low a value for the fitted exponent. As 
with the data of Bastea fc Lebowitz (1997 ) and Laradji et al. (1996 ), the L(T) plot is 
lacking the initial diffusive plateau, an absence that we found was invariably associated 
with strong residual diffusion within LB. 



B.4. \Jury, Bladon, Krishna & Cotes (1999\\ ) 



Jury et al. (lyyyb) carried out a series of simulations of a symmetric, binary fluid mix- 
ture using the DPD (dissipative particle dynamics) method with 10 6 particles and a 
deep quench. (The DPD algorithm combines soft interparticle repulsions with pairwise 
damping of interparticle velocities and pairwise random forces.) In terms of the range 
of domain scales that can be probed, as a multiple of the interfacial width, these are 
roughly equivalent to our 128 3 runs using LB. These authors found that each data set 
was well fitted by a linear scaling, L = v(T — Ti nt ), but with a systematic increase of 
b± = vj (Lq/Tq) upon moving from upper right to lower left in the universal scaling plot, 
see figure [l9| (inset). These seven data sets have a range of L values, 0.29 ^ L ^ 0.013, 
which places them all firmly in the crossover region found in the LB results, between 
Run029 (L = 0.95) and Run030 (L = 0.01). 

Jury et al. suggested two alternative interpretations of their own data and that of 



Laradji et al. (1996) and Bastea & Lebowitz (1997), to explain the observed linear scaling 
within each run, but lack of consistency in the prefactor, b%. The first was a possible 
nonuniversality arising from the physics of pinchoff. The second was that all data sets 
formed part of an extremely broad crossover region. This explains the observed trend 
bi r~j i~ ' 2 , but not the linear scaling within each run. However, Jury et al. noted the 



possibility of 'dangerous' finite size effects within their data, which was truncated at 
L m ax = A/2. Our LB results support the latter explanation for the DPD data, with a 
separate or additional reason (residual diffusion) for rejecting the other datasets. Unlike 
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these, all the data sets of Jury et al. do lie very close to the LB results, see figure 19 
(inset). Moreover, the DPD show a diffusive plateau prior to the onset of the apparent 
linear regime, consistent with low residual diffusion levels. 

Since the DPD simulation method is very different from LB, the correspondence of the 
two sets of results lends broad support to the idea of a universal scaling, although the 
fact that each DPD run is best fit by a locally linear growth law is not consistent with 
this. Based on our LB results, we are inclined attribute the latter to finite size effects, 
which would certainly be large enough to spoil our nearest equivalent runs (128 3 with 
Lmax = 64). However, we cannot rule out other explanations as offered by Jury et al 



Larger DPD runs, and/or a series of LB runs, giving strongly overlapping datasets in a 
narrow window of the l(t) like the DPD results (figure 19, inset), might shed light on 
this. 
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